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Steady,  transonic  flow  over  a  projectile  was  numerically 
investigated  to  gain  insight  and  understanding  of  the 
thin-layer  Navier-Stokes  approximation  and  to  assess  the 


forces  acting  on  a  projectile  at  transonic  speeds.   The 
turbulence  model  used  in  the  computation  is  an  algebraic 
two- layer  eddy  viscosity  turbulence  model  proposed  by 
Baldwin  and  Lomax.  The  projectile  model  considered  was  a 
secant-ogive-cylinder-boattail  configuration  for  which 
accurate  surface  pressure  measurements  are  available  to 
assess  the  numerical  results.  The  sensitivity  of  the 


relative  importance  of  viscous  effects  on  the  aerodynamic 


numerical  scheme  to  the  grid  resolution  and  added  numerical 
dissipation  terms  was  investigated.  The  surface  shear  stress 
distribution  has  also  been  computed  to  aid  in  identifying 
regions  of  the  reversed  flow. 

The  computed  results,  based  on  a  grid  system  generated  by 
a  set  of  hyperbolic  partial  differential  equations,  were 
verified  by  comparing  with  experimental  measurements.   Good 
agreement  between  the  numerical  results  and  measured  data 
indicates  that  the  numerical  scheme  can  predict  the 
essential  features  of  the  flow  field.   The  calculated 
results  obtained  at  zero  angle  of  attack  and  Mach  number 
0.91  show  that  the  viscous  drag  is  as  important  as  the 
pressure  drag  (no  base  drag) .  However,  the  relative 
importance  of  viscous  drag  decreases  with  increasing  Mach 
number  and  angle  of  attack.  At  10  and  45  degrees  angle  of 
attack,  the  viscous  forces  had  a  negligible  effect  on  both 
the  lift  force  and  the  pitching  moment. 
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CHAPTER  I 
INTRODUCTION 


An  accurate  prediction  of  the  aerodynamic  forces  is 
important  in  the  design  of  flight  vehicles.  For  complex  flow 
fields,  an  accurate  simulation  of  the  flow  field  is  crucial 
for  the  successful  prediction  of  such  forces.  Considerable 
experimental  efforts  have  been  undertaken  to  enhance  our 
understanding  of  these  complex  flows  and,  recently,  with  the 
increase  in  computational  capabilities,  there  has  been  a 
major  effort  to  numerically  simulate  such  flows.  The  most 
challenging  features  of  transonic  flow  to  be  numerically 
simulated  are  the  viscous/inviscid  and  shock/boundary  layer 
interactions.  There  are  two  basic  approaches  used  to 
numerically  simulate  these  flows.   One  method  is  a  patching 
technique  in  which  the  solution  is  obtained  from  inviscid 
equations  and  boundary  layer  equations  separately  and  then 
patched  together  in  some  iterative  process.   Another 
approach  is  to  solve  a  suitable  subset  of  the  Navier-Stokes 
equations  for  the  entire  flow  domain.  The  primary  advantage 
of  the  patching  method  is  their  computational  efficiency. 
However,  it  cannot  solve  the  separated  flow.  The  latter 
approach  has  the  advantage  that  it  provides  a  continuous 
treatment  of  the  solution  through  the  interaction  regions. 
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In  the  transonic  flow  regime,  critical  aerodynamic 
behavior  has  been  observed  which  can  be  attributed  to  the 
complex  shock  pattern  and  local  flow  expansions  that  exist 
over  the  projectile  at  these  speeds.  Shadowgraphs  which  can 
be  found  in  reference  7,  show  the  shock  pattern  existing  on 
a  typical  projectile  configuration  for  transonic  flow.  The 
formation  of  shock  waves  imbedded  in  the  flow  field  produces 
severe  changes  in  the  aerodynamic  coefficients  such  as  the 
drag  force  and  the  pitching  moment.  For  example,  the 
nondimensional  drag  force  for  a  typical  projectile  shape  has 
been  found  to  change  by  as  much  as  100%  through  a  Mach 
number  range  of  0.95  to  0.97  [9].  Similar  results  have  also 
been  evident  in  this  study.  For  such  large  changes  in  the 
^magnitude  of  the  aerodynamic  forces,  it  is  essential  to 
understand  and  be  capable  of  computing  the  features  of  the 
flow  field  that  contribute  to  this  phenomenon.  In  the 
transonic  flow  regime,  the  flow  field  is  characterized  by 
strong  viscous/inviscid,  shock/boundary  layer  interactions. 
It  is,  therefore,  advantageous  to  use  the  approach  which 
solves  the  suitable  subset  of  Navier-Stokes  equations  since 
this  method  considers  these  interactions  in  a  fully  coupled 
manner. 

A  three  dimensional  code  was  developed  in  1978  by  Pulliam 
and  Steger  at  NASA  Ames  to  solve  the  thin-layer 
Navier-Stokes  equations  for  an  ideal  gas  in  a  transformed 
boundary- fitted  coordinate  system  [14].  In  this  code,  which 
is  applicable  to  high  speed  compressible  flow  problems,  the 
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governing  equations  are  approximated  by  the  spatially 
factored  implicit  scheme  of  Beam  and  Warming  [3,4].  To 
control  the  numerical  stability  of  the  solution  algorithm, 
second  order  implicit  and  fourth  order  explicit  dissipation 
terms  are  added.  Consequently  it  results  in  a  block 
tridiagonal  system  of  equations  to  be  solved  at  each  time 
level.  The  turbulence  closure  model  implemented  in  the  code 
is  an  algebraic  two-layer  eddy  viscosity  model  of  Baldwin 
and  Lomax  [2].   The  unsteady  thin-layer  Navier-Stokes  code 
has  an  option  for  calculating  inviscid  flow  problems  and  a 
steady  flow  solution  can  be  obtained  as  the  converged 
solution  of  a  time  asymptotic  flow  problem.  This  code  has 
been  widely  applied  to  various  flow  problems. 


The  three-dimensional  thin-layer  Navier-Stokes  code  has 
also  been  simplified  for  axisymmetric  flow  cases.  Both  codes 
have  been  investigated  to  some  extent  by  the  U.S.  Army 
Ballistic  Research  Laboratory  on  transonic  projectile 
aerodynamic  problems  [9,11].   The  planar  grid  network 
provided  to  a  Navier-Stokes  code  is  obtained  from  a  grid 
generation  code  GRIDGEN  which  can  give  either  an  elliptic 
grid  or  a  hyperbolic  grid  [10].  Moreover,  a  three 
dimensional  grid  system  is  formed  by  a  sequence  of  planar 
grids  around  the  axis  of  a  projectile.  For  a 

secant-ogive-cylinder-boattail  projectile  with  sting  at  zero 
angle  of  attack,  the  published  result  showed  that  the 
computed  surface  pressure  coefficient  over  the  secant-ogive 
portion  and  the  boattail  portion  of  the  projectile  agrees 
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rather  well  with  the  measured  data  but  the  agreement  on  the 
cylinder  portion  of  shock/boundary  layer  interaction  region 
is  not  very  satisfactory  for  some  cases  considered.  For  the 
projectile  model  at  2  degrees  angle  of  attack  the 
Cp-distribution  computed  agrees  qualitatively  with  the 
measured  data  but  quantitatively  the  agreement  over  the 
cylinder  portion  and  boattail  portion  is  not  satisfactory 
[18].   The  published  results  seem  to  indicate  that  the 
thin- layer  Navier-Stokes  codes  can  give  acceptably  accurate 
surface  pressure  for  the  complex  transonic  projectile 
aerodynamic  problem  if  a  good  adaptive  grid  system  is 
provided.  However,  more  precise  causes  for  the 
unsatisfactory  results  reported  are  yet  to  be  investigated; 
moreover,  no  result  on  the  skin  friction  coefficient  Cf  has 
been  reported  and  discussed. 

The  GRIDGEN  code  obtained  from  the  U.S.  Army  has  been 
investigated  and  modified  [6].  The  modified  hyperbolic  grid 
is  generated  using  predetermined  adaptive  boundary  points. 
In  order  to  gain  the  smoothness  property,  the  hyperbolic 
grid  obtained  from  GRIDGEN  also  has  been  modified  to  bend 
the  grid  lines  in  the  nose  region  of  the  projectile,  which 
has  been  found  to  be  important  in  the  convergence  process 
and  for  the  solution  accuracy.  The  modified  hyperbolic  grid 
generation  code  is  efficient,  requiring  less  computer  time 
than  elliptic  grid  generation  techniques;  moreover,  it  can 
provide  a  good  grid  network  and  is  therefore  used  in  this 
study. 
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An  accurate  prediction  of  a  complex  viscous  flow  field  via 
thin-layer  Navier-Stokes  codes  depends  upon  many  factors. 
These  include  grid  resolution,  turbulence  models,  artificial 
dissipation  terms,  flow  characteristics,  etc.  For  a 
successful  application  and/or  modification  of  a 
Navier-Stokes  code,  one  must  be  familiar  with  certain 
relevant  features  of  the  code.  This  motivates  the  present 
investigation  of  the  computer  codes.  The  main  objectives  of 
this  study  are  to  further  advance  our  understanding  and  gain 
in-depth  insight  on  the  application  of  thin-layer 
Navier-Stokes  codes  implemented  with  the  Baldwin-Lomax 
turbulence  model  for  accurate  computation  of  transonic 
projectile  aerodynamics. 


CHAPTER  II 
GOVERNING  EQUATIONS 


Navier- Stokes  Equations 

In  three  dimensions,  the  unsteady  compressible 
Navier-Stokes  equations  referenced  to  Cartesian  coordinates 
without  body  forces  and  external  heat  addition  can  be 
written  in  the  following  vector  form: 


3q    3E    3F    3G    3#v    3FV    3G 

+  —  +  +  —  =  +  +  (l) 

3t    3x    3y    3z    3x     3y     3z 


in  which  q   is  an  unknown  vector,  E,    F,    G   are  the  inviscid 
flux  vectors  and  Ev    ,  Fv  ,  Gv  are  the  viscous  flux  vectors. 
The  explicit  expressions  for  these  vectors  can  be  written  as 


m 

q   =    [    p  pu  pv  pw  E  ] 

2  T 

E  =    [    pu  pu  +p  puv  puw  (E+p)u  ] 

F  =    [    pv  puv  pv  +p  pvw  (E+p)v  1 

G  =    [    pw  puw  pvw  pw  +p  (E+p)w  J  (2) 

£v  =  I  °  Txx  Txy  *xz  ^x  '* 
Fv  "  t  °  Tyx  Tyy  Tyz  Py  ]T 
Gv  =  [  °  *zx  Tzy  Tzz  *z  1T 
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Here,  p  is  the  density;  u,  v  and  w  are  the  velocity 
components  in  x,  y  and  z-directions,  respectively;  p  is  the 
pressure;  and  E  is  the  total  energy  per  unit  volume  which  is 
related  to  the  internal  energy  e  and  the  kinetic  energy  by 
the  formula 


1 

E  =  p  [e  +  -  (u^  +  vz  +  wz)]  (3) 

2 


For  a  Newtonian  fluid  with  zero  bulk  viscosity,  the 
relationships  between  stresses  and  velocity  gradients  are 


2  3u    3v    3w 
xxx=-y(2  ) 

3  3x    3y    3z 

2  3v    3w    3u 

Vy  =  _li(2 ) 

3  3y    3z    3x 

2  3w    3u    3v 
xzz=-U(2  ) 

3  3z    3x    3y  (4) 

3u    3v 

Txy  =  ii  (  —  -—)  =  x 

3y   3x      J 

3v    3w 

Tyz  =  y  ( )  =  Tzy 

*  3z    3y        y 

3w    3u 

Tzx  =  "  (  —  -  —  )  =  Txz 
3x    3z 


and  the  3 ' s  in  equations  (2)  are  related  to  viscous 
dissipation  and  heat  conduction  by 
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3T 

Px   =  UTxx   +   VTxy  +   WTxz    +    k  — 

3x 

3T 
eY  =  UTyx  +  VTyy  +  WTyz  +  k  ~  (5) 

3T 

Pz  =  UTzx  +  VTzy  +  wtzz  +  k  — 

3z 


where  y  and  k  are  the  coefficients  of  dynamic  viscosity  and 
thermal  conductivity,  respectively,  and  T  is  the 
temperature. 

The  complete  set  of  time-dependent  Navier-Stokes  equations 
(1)  includes  one  continuity  equation,  three  momentum 
equations  and  one  energy  equation,  but  there  are  seven 
unknowns  p,  u,  v,  w,  E,  p  and  T.   In  order  to  close  the 
system  of  equations,  the  fluid  is  assumed  to  be  a  perfect 
gas.   The  equation  of  state  and  other  relations  for  a 
perfect  gas  are 


p  =  PRT 


e  =  cvT   h  =  c  T    %    = 


cp  (6) 


cv 


where  R  is  the  gas  constant,  h  is  enthalpy,  and  c   and  c 

■        p      v 

are  the  specific  heat  at  constant  pressure  and  volume, 
respectively.   JT  is  the  ratio  of  specific  heats,  which  for 
air  at  standard  conditions  is  1.4.  Now,  by  using  equations 
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(3)  and  (6),  one  can  obtain  two  additional  equations 
relating  p  and  T  to  the  dependent  variables  q 


p  =  (Z    -  1)[E  -  -  P  (u2  +  v2  +  w2)] 

2 


(Y-l)   El 

T  =  [ (u2  +  v2  +  w2)] 

R     p    2 


(7) 


Hence,  equations  (1)  and  (7)  are  the  governing  equations  for 
a  class  of  unsteady  compressible  flow  problems. 

The  number  of  parameters  for  flow  problems  of  the  same 
class  can  be  reduced  by  introducing  dimensionless  variables. 
Once  the  equations  are  put  into  nondimensional  form,  the 
solution  can  easily  be  used  to  correlate  data  for  succinct 
presentation  using  the  minimum  possible  amount  of 
experimental  testing.  The  dimensionless  quantities  used  here 
are  obtained  by  choosing  characteristic  scales  in  the 
following  manner: 


* 

X 

* 

y 

y 
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z 

z 
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t 
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L 

L 

L 
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(8) 
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where  the  nondimensional  variables  are  denoted  by  an 
asterisk  *  and  the  freestream  conditions  by  ~ .  For  instance, 
VM  is  a  characteristic  velocity  of  the  fluid  at  freestream 
and  L  is  the  characteristic  length  scale. 

The  substitution  of  the  dimensionless  quantities,  equation 
(8)  into  the  governing  equations  (1)  and  (7),  yields 


3q*    3E*    3F*    3G*  3B  *    3F  *    3G  * 

— *  +  — *  +  — *  +  —  =  Re      i—ir  +  —T-  +  — T")  (9> 

St     3x     3y     3z  3x      3y      3z 


and 


p*  -  (T  -  1)[E*  -  -  P*  (u*2  +  v*2  +  w*2)] 

2 

(10) 

E     1 

T*  =  1(1-1)  Mj    [—  -  -  (u*2  +  v*2  +  w*2)] 

P     2 


where  the  Reynolds  number  Re  and  the  Mach  number  at 
freestream  are  defined  as 


OO        00  oo 

Re   =  Mm   =       (11) 

y_  /2TRT 


The  unknown  vector  q      and  flux  vectors  E* ,  F* ,    G* ,    E   *, 
Fv    >    Gv     have  the  same  form  of  equation  (2)  as  well  as  the 
stresses  in  equations  (4)  except  that  all  the  quantities  are 
nondimensional  quantities.  The  components  of  &  become 


-li- 
ft * 

„*******  !  J 

PV        =    U     TvV        +     V     T„„        +     W     t„_         + 


x  xx  'xy  "    lxz       "  2        * 

(J-l)Pr  M/    3x 

*  * 
„*            ******                       V                  3T 

Py        =    U     X  +VT  +WT  f     -_  (12) 

1  xx  *  (y-l)Pr  M/    3y 

*  * 
_*******                   v               9T 

Bz      =UIzx      +   v   Tzy      +WTzz      +   5  ~ 

1  (I-l)Pr  M  Z  3z 


Here,  the  Prandtl  number  Pr  is 


c_  u 

Pr  =  — - 


(13) 


The  typical  values  of  the  Prandtl  number  for  air  at  standard 
conditions  are  0.9  and  0.72  for  turbulent  and  laminar  flow, 
respectively.  In  order  to  simplify  the  notation,  the 
asterisks  are  removed  from  the  nondimensional  equations 
hereafter  except  where  noted. 

Transformed  Navier-Stokes  Equations 

To  enhance  numerical  simplicity,  accuracy  and  efficiency, 
a  generalized  boundary-fitted  curvilinear  coordinate  system 
is  used.   The  boundary-fitted  coordinate  transformation  maps 
all  boundary  surfaces  in  the  physical  domain  onto  coordinate 
surfaces  of  the  computational  domain.  The  boundary- fitted 
curvilinear  coordinates  not  only  facilitate  implementation 
of  boundary  conditions  accurately,  but  also  allow  all 
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computations  to  be  performed  on  a  fixed  square  grid  in  the 
computational  domain.  The  numerical  solution  technique, 
thus,  becomes  more  versatile. 

Consider  a  generalized  transformation  of  the  form 

£  =  £  (  x,y,z,t  ) 

n  =  ti  (  x,y,z,  t  ) 

C  =  C  (  x,y,z,t  ) 
x  =  t 

where  (£,ti,C,t)  are  the  coordinates  of  the  computational 
domain  and  (x,y,z,t)  are  those  of  physical  domain.  It  is 
known,  e.g.  reference  1,  that  the  Jacobian  of  the 
transformation  J~  ,  which  represent  the  volume  of  the  grid 
cells  in  physical  space,  is  given  by 


J_1  =  x£ynzc  +  Vc2*  +  xsVti  "  xSyCzn  "  xcynzS    <14> 


and  the  metrics  in  the  computational  domain  can  be  related 
to  those  in  the  physical  domain  by  using  a  chain  rule 
expansion  of  x^ ,  y^ ,  etc.,  and  then  solving  for  £x/  £  , 
etc . ,  to  give 
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*x  =  J  (  y,zc  -  yczn 

KY  =   J    (    z^    -    z^ 

*z  =  J  (  *,yc  -  x^t, 


=   J    (    yrz 


n,,  = 


czS  "  ^zc 

J    (    zcxc    -    z^xc 


=   J    (    x    y,    -    x,y 


Z*Z 


'KH 


^x   = 


C.    = 


u  = 


J   (   nzn   "  V« 

J   (   ZK\   "  zr,xC 
J  {  xSyr,  "  xi,y5 


"Mx    "    My    "    M: 


T1+.     =     "XJ 


«t  = 


t-x    "    y^y    -    ZTT1Z 
"Mx   "   yT5y  -   zt^z 


Note  that  the  subscript  in  the  metrics  denotes 
differentiation.  The  metric  identities  of  the 
transformation,  e.g.  reference  1  or  20,  are 


3    1     3    Kt  a    nt     8    Ct 

—  (-)+—(  —  )+—(  —  )  +  —  (  —  ) 

3t    J      3$    J       8ti    J       3C    J 

9     Cx      3     ^x      3     Sx 

—  (  —  )+—(  —  )+—  (  -^  ) 
3?    J       3ti    J       3£    J 

3  ,  *y\   a    "y    3    cy 

—  (  —  )+  —  (  -*■  )  +  —  (  -i  ) 


3?    J 

3    5, 


3n   J 

3     ti, 


3C    J 
3     C. 


3£    J       3ti    J       3C    J 


=  0 


=  0 


0 


=  0 


(15) 
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By  applying  the  generalized  transformation  to  the 
compressible  Navier-Stokes  equations  (9)  in  vector  form, 
dividing  by  the  Jacobian  and  then  rearranging  into 
conservation  law  form  with  the  use  of  metric  identities, 
equation  (15),  the  following  form  of  the  transformed 
Navier-Stokes  equations  is  obtained: 


3q    3E    3F 

—  +  —  +  +  — 

3t    3£    3T)    3£ 


3G       .    3E 
-1      v 

=  Re  1    (  — - 


3Fv    9Gv 

+      +  )     (16) 


3?     3n 


3C 


in  which 


-1-1  T 

q  -  J  q  =   J   [  p  pu  Pv  pw  E  ]x 


_  t-1 


3£     3C     3?     3£ 
E  =  J  x  (  — q   +   — E   +   — F   +  — G    ) 
3t-    3x     3y     3z 


-    T-1 


-J    [  PU  PUu+£xp  pUv+Cyp  pUw+CzP  (E+p)U-£tp  J 
3ti     3t|     3  n     3t) 


F  =  J-l  , 


r-l 


— q   +   — E   +   — F   +   — G    ) 
3t     3x     3y     3z 


-  J    [  pV  pVu+nxp  pVv+tiyp  pVw+nzp  (E+p)V-ntp  ]T 

,    U     3£     3C     3£ 
S  =  J    (  — q   +  — E  +   — F  +   — G    ) 
3t     3x     3y     3z 

=  J"   [  PW  PWu+Cxp  pWv+Cyp  pWw+CzP  (E+p)W-Ctp  ]T 

0 

?xTxx  +  Vyx  +  ?zTzx 


?xTxz  +  *yTyz  +  CzT 
*xPx  +    *y*y  +    5ZPZ 


+   £   T       +   £   T 

xy   *y  yy   ^zTzy 


V 
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nxTxx  +  Vyx  +  'Vzx 


F„  = 


J 


G„  =  - 


+  Tl  T     +  n  T 

x  xy    'y  yy    'z^zy 


Vxz  +  Vyz  +  'Vzz 

Vx  +  Vy  +  .Vz 


f 


sx  xx  ^ylyx    ^z  zx 

**xlxy  ^y  yy    4Tzy 

*x  xz  y  yz   ^z  zz 

«x*x  +  Sey  +  cz*z 


y 


The  stresses  t  and  temperature  gradients,  in  the  p  terms 
contain  the  metrics  of  the  transformation;  e.g., 


3u  3v     3u    3u    3u    3v    3v    3v 
Txy=lJ< )=w(— 5  y+— n  +— C Cx nx ^x)=txv 

3y  3x     3?  Y  3t,  Y  3?  Y  3^    3n    3C       Y 


etc.  U,  V  and  W  are  the  contravariant  velocity  components 
along  the  three  curvilinear  coordinates  £,  ti  and  c, 
respectively.  They  are  defined  in  terms  of  the  Cartesian 
velocity  components  u,  v  and  w  as 


U   =    *t   +    U^x    +    v*y   +    w^z 
V  =  nt  +  utix  +  vti      +  wtiz 

W  =    Ct   +   ucx   +    vcv   +    wcz 


(17) 


-  16  - 

Thin  Layer  Navier-Stokes  Equations 

In  general,  a  very  large  amount  of  computer  time  and 
extensive  storage  are  needed  for  solving  a  complex 
aerodynamic  problem  governed  by  the  complete  Navier-Stokes 
equations  (16).  At  high  Reynolds  numbers,  however,  the  flow 
around  the  body  can  be  considered  to  be  mostly  inviscid, 
except  in  a  thin-layer  region  close  to  the  body  surface 
where  the  viscous  effects  become  important.   With  this  in 
mind,  the  thin-layer  approximation  is  used  to  simplify  the 
Navier-Stokes  equations  in  generalized  coordinates. 
The  thin-layer  approximation  to  the  Navier-Stokes 
-  equations  is  based  on  an  assumption  that  the  diffusion 
processes  along  the  body  surface  are  negligibly  small  in 
comparison  with  the  diffusion  process  in  the  direction 
normal  to  the  body  surface.  This  is  generally  true  for  high 
Reynolds  number  flow.   As  a  consequence  of  this  assumption, 
all  viscous  terms  containing  derivatives  parallel  to  the 
body  surface  are  dropped  since  they  are  substantially 
smaller  than  viscous  terms  containing  derivatives  normal  to 
the  surface.   Assume  that  a  body  surface  is  mapped  into  a 
£ti -plane  in  the  transformed  space.  Then  £  and  r\    derivatives 
of  the  viscous  terms  can  be  neglected.  Hence,  the  thin- layer 
Navier-Stokes  equations  in  curvilinear  coordinates  are 
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3S 


-1 


3q    3E    3F    3G 

—  +  +  +  =  Re""1- 

3t    3£    3  ti    3£         3C 


(18) 


where 


f 


S  =  J"1 


mlUC  +  m2«x 
mlvC  +  m2^y 
mlwC  +  m2^z 


vmlm3    +    m2(Cxu   +    SV   +    CzW)/ 


ml    =   *<<x      +    V    +    Cz    > 

m2  =  y(Cxur  +  cvvr  +  Cwwr)/3 


'x"c 
m3  =  (u2  +  v2  +  w2)f/2  +  Pr"1(y-l)"1(c2) 


•y  c 


and  c=  /2TRT  is  the  local  speed  of  sound.  Since  the 

thin- layer  Navier-Stokes  equations  (18)  retain  the  normal 

momentum  equation,  the  pressure  can  vary  through  the  viscous 

layer.   Hence,  the  thin- layer  model  is  more  general  than  the 

boundary  layer  equations  and  is  capable  of  predicting  flow 

separation. 

Axi symmetric  Thin-Laver  Navier-Stokes  Equations 


For  an  axisymmetric  flow  in  which  the  flow  variables  are 
invariant  in  the  n-direction,  the  three-dimensional 
thin- layer  Navier-Stokes  equations  can  be  further  simplified 
by  imposing  the  following  restrictions:  (i)  the  body 
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geometry  is  axi symmetric;  (ii)  the  flow  variables  and  the 
contravariant  velocities  do  not  vary  in  the  ri-direction. 
Hence,  the  3F/3ti  term  of  equation  (18)  can  be  reduced  to  a 
source  term  H.   Based  on  the  relationships  between  the 
coordinate  systems,  the  Cartesian  coordinates  (x,y,z),  the 
cylindrical  coordinates  (x,<fi,R)    and  the  curvilinear 
coordinates  (£,t),c)  shown  in  figure  1,  one  has 

x  =  x(  s,c,t  ) 

y  =  R(  £,C,t  )  sintf 

Z  =  R(   $,S,T   )   COS0 

By  evaluating  the  metric  terms  with  the  use  of  the  above 
relations  and  substituting  them  into  transformed  thin-layer 
Navier-Stokes  equations  (18),  the  resulting  unsteady 
axi symmetric  thin- layer  Navier-Stokes  equations  become 


3q    3E    3G  3S 

—  +  —  +  —  +H=  Re-1  —  (19) 

3x    3£    3C  3C 


With  the  metrics  and  flux  vectors  evaluated  at  plane  <t>    =   0° , 
the  source  vector  H  [9,11]  can  then  be  expressed  as 
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Projectile-type  body 


x 


-  constant  plane 


0  =  constant  plane 


Figure  1.  Axisynunetric  body  and  coordinate  system 


H 


-  20  - 

0 
0 
J'1  (  PV[R5(U-5t)  +  R5(W-«t)] 
-PVR  (V-nt)  -  p/R 
0 


-  t-1 


which  has  replaced  3F/3 n  of  equation  (18). 

Equation  (19)  contains  only  two  spatial  derivatives  but 
does  retain  all  three  momentum  equations  by  allowing 
non-zero  velocity  components  in  the  invariant  Ti-direction. 
Thus,  flow  over  axisymmetric  spinning  projectiles  can  be 
solved  with  equation  (19). 


CHAPTER  III 
NUMERICAL  ALGORITHM 


There  have  been  various  numerical  algorithms  developed 
for  approximating  the  Navier-Stokes  equations.   For  the 
thin-layer  Navier-Stokes  codes  obtained  from  the  U.S.  Army 
Ballistic  Research  Laboratory  and  used  in  this  study, 
equations  (18)  and  (19)  are  approximated  by  the  Beam  and 
Warming  scheme  [3,4,21],  which  is  an  implicit  approximate 
factorization  finite  difference  scheme.   This  effective 
finite-difference  scheme  is  described  and  discussed  in  the 
following.  Note  that  for  a  fixed  grid  system,  the  t  equals 
to  t. 


Time-Differencing 

Consider  three  different  temporal  differencing 
approximations  given  by 

(1)  Backward  Euler  implicit  formula 


qn+1  -  q11  =  At  (  — )n+1  +  0[(At)2] 

at 


21  - 
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(2)  Trapezoidal  formula 


.  At    3q        3q 

qn+l  .  qn  =  _  [(_}n  +  (_)n+l]  +  0[(At)3] 

2    at     at 


(3)  Three-point  backward  formula 


3qn+1  -  4qn  +  qn_1  =  2At  (  —  )n+1  +  0[(At)3] 

at 


By  defining  Aqn  =  q11    -  qn ,    the-  above  formulas  can  be 
combined  into  a  single  generalized  time-differencing  formula 
in  delta  form 


aq  _..  aq 

rAt  =  (1  +  s)Aqn  -  sAq11"1  -  At  ( 1  -  r) 

at  at 

+  0[ (r-s-l/2)(At)2  +  (At)3]  (20) 


With  the  appropriate  choice  of  the  parameters  r  and  s, 
equation  (20)  reproduces  many  familiar  two-  and  three- level 
explicit  and  implicit  schemes  which  are  listed  in  the 
following  table. 
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Table  1.  Partial  list  of  schemes  contained  in  equation  (20). 

r     s        Scheme  Truncation  error 


0     0  Euler,  explicit  0[At] 

0  -1/2  Leapfrog,  explicit  0[(At)2] 

1  0  Euler,  implicit  0[At] 
1/2     0  Trapezoidal,  implicit  0[(At)2] 

1      1/2  3-point-backward,  implicit  0[ (At)  ] 


After  substituting  equation  (18)  into  equation  (20),  we  have 


rAt  3E    3F    3G       , 3S    ,     s 

Aq11  +  (—  +  —  + Re"1— )n+1  =  Aq11"1 

1  +  s  3£    3l|    3£        ac        i  +  s 

1-r  3E    3F  3G         3S 

-At  ( —  +  —  + Re"1 — )n  (21) 

1  +  s  3$    3t)  3£        3£ 


which  is  either  first  order  or  second  order  in  time, 
depending  on  the  choice  of  r  and  s  as  indicated  in  the  above 
table . 

The  flux  vectors  E,  F  and  G  at  the  advanced  n+1  time  step 
are  nonlinear  functions  of  the  unknown  vector  q.  It  is 
possible  to  remove  the  nonlinear  nature  and  still  maintain 
the  order  of  accuracy  of  the  finite  difference  approximation 
by  using  Taylor  series  expansion  about  q11 
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En+l  =  En  +  AnAqn  +  0[(At)^j     where   A  =  — 

3q 

3  F 
Fn+1  _  pn  +  BnAqn  +  0[(At)2]     where   B  =  — 

3q 

3G 
Gn+l  _  Gn  +  cnAqn  +  0[(At)2]     where   C  =  — 

3q 


Matrices  A,  B  and  C  are  so  called  Jacobian  matrices.   The 
viscous  flux  term  can  also  be  linearized  by  using  a  method 
suggested  by  Steger  [16].  In  order  to  apply  this 
linearization  method,  the  coefficients  of  viscosity  ]i    and 
thermal  conductivity  k  are  assumed  to  be  locally  independent 
of  q.  As  a  result  of  this  assumption,  elements  of  S  have  the 
general  form 


_!     9 

si   =    J         ai   —  (»*> 

3C 


where  a-  is  independent  of  q  =  J~  q   given  in  equation  (16), 
and  B-  is  a  function  of  q.   These  elements  are  linearized  in 
the  following  manner  by  Taylor  series  expansion 


s;n+1  =  s/1  +  J"1  a/1  —  [I    (—V  Aqmn]  ♦  0[(At)2] 

3C  m      3qffl 
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These  relations  can  be  written  in  a  vector  form  as 


;n+l  =  sn  +  j-lMnJAqn 


which  has  a  truncation  error  of  second  order  in  time.   The 
explicit  expression  of  Jacobian  matrices  A,  B,  C  and  M  is 
given  in  the  Appendix.  Equation  (21)  now  becomes 


rAt  3A    3B    3C         3 

[  I  +  ( +  + Re"1 (J_1MJ)]nAqn 

1+S  3£    3n    3£         3£ 

s  At  3E  3F    3G        3S 

Aq11"1 ( +  + Re"1 — )n      (22) 


1+s  1+s  3?    3n    3£         3£ 


where  the  matrix  I  is  the  identity  matrix. 


Approximate  Factorization 

To  increase  computational  efficiency,  approximate 
factorization  is  employed  for  the  solution  algorithm  of 
equation  (22),  which  reduces  the  solution  process  to  three 
one-dimensional  problems  at  a  given  time  level.  The  solution 
process  for  Aq11  then  involves  the  inversion  of  three 
matrices  with  smaller  bandwidth,  for  which  there  exist 
efficient  solution  algorithms.  Consider  the  following 
factorized  equation  approximating  equation  (22). 
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rAt  3A  rAt  3B 

[  1+ ]n[  1  + ]n 

1+s  35  1+s  3ti 

rAt  3C         3 

[  I  + Re"1  — (J_1MJ)  ]n  Aq11  =  RHS(22)   (23) 

1+s  3£         3C 


Here  RHS(22)  is  used  to  indicate  the  right-hand  side  of 
equation  (22).   The  error  introduced  by  the  factorization 
procedure,  equation  (23),  is  the  leading  third  order  term 


At3  3A  3B  3B  3C  3C  3A   At4  3A  3B  3C    3qn 

[ ( + + )  + ]n  +0[  (At)5] 

4   3^  3  n  3n  3£  3C  35    8   35  3ti  3C    3t 


which  does  not  affect  the  second  order  accuracy  of  the 
difference  approximation. 


Numerical  Dissipation 

For  a  complex  aerodynamic  problem,  a  direct  application 
of  equation  (23)  often  experiences  convergence  problems. 
Hence,  dissipation  terms  are  added  to  the  equations  to  damp 
out  high  frequency  oscillations  and  to  control  numerical 
instabilities.  Several  different  numerical  dissipation  terms 
have  been  investigated,  e.g.  reference  13.  Currently,  the 
numerical  dissipations  proposed  by  Beam  and  Warming  [3,4] 
and  Steger  [16]  are  employed  in  the  thin-layer  Navier-Stokes 
codes.   The  second  order  implicit  dissipation  operators 


-  27  - 

DH  =  ~£i  J~\V 

DH   =  "El   J'\V  <24> 

DU  =  "ei  J"lvcV 


are,  respectively,  added  to  the  left  hand  side  implicit 
one-dimensional  £ ,  n  and  c  operators  of  equation  (23)  while 
an  explicit  fourth  order  dissipation  operator 


DE  =  -tE    J_1[(VCAC)2  +  (V/^)2  +  (V^)2]J       (25) 


is  added  to  right  hand  side  of  equation  (23).  The  symbols  A 
and  V  are  the  standard  forward  and  backward  difference 
operators,  that  is, 


1 
h 


™i  =  -  (fi  -  fi-0  +  °(h) 

h 


Finite  Difference  Equations 

The  spatial  derivatives  appearing  in  equation  (23)  have 
to  be  represented  by  finite  difference  quotients.  In  the 
computer  codes,  central  difference  approximations  are  used. 
Thus,  the  equation  (23)  is  rewritten  as 
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5       H       C 


L?    L      L_    Aq11   =    Rn  (26) 


where 


rAt 

Lr    =    I    +   6?An   +    DTr 

*  l  +  s      *  ift 

rAt 

Lr,    E    x    +    6„Bn    +    DT„ 

l+s 

rAt  rAt 

L      =    I    +    6_Cn    Re_16     (J_1MnJ)    +    DT 

g  l+s      g  l+s  g  * 

s  At 

Rn    =    Aq11"1 (6rE    +    6    F    +    5rG    -    Re_16_S)n    +    Dvqn 

l+s  l+s 


At  this  point,  equation  (26)  can  be  first  or  second  order 
accurate  in  time  depending  on  the  choice  of  parameters  r  and 
s.   For  the  finite  difference  approximation  of  spatial 
derivative  terms,  a  second  order  central  differencing 


*fi   =~    <f;+7  "  fi-0    +    °<h2) 
2h 


is  employed  on  the  implicit  part  of  fi^A,  6  B  and  6_C  for 
maintaining  the  block  tridiagonal  matrix  structure,  while  a 
fourth  order  central  difference  approximation 


6f*  =  7^  ("f**  +  Sf^  '  8f^-'  +  £*-*}  +  0(h4) 

12h 
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is  used  on  the  convection  terms  6>.E,  6  F  and  6„G  of  the 
RHS(26)  to  keep  the  convection  truncation  error  from 
exceeding  the  magnitude  of  the  truncation  error  of  the 
viscous  terms  6  S  which  is  approximated  by  second  order 
central  difference.   In  this  study,  the  Euler  implicit 
two-level  scheme  is  used;  i.e.,  r=l  and  s=0.  Therefore,  the 
solution  algorithm  is  first  order  accurate  in  time  and 
second  order  accurate  in  space. 

The  unknown  vector  q  at  n+1  time  level  governed  by 
equation  (26)  can  now  be  computed  effectively  as  follows: 


**     n 
L?  Aq    =  Rn 

•k  *  -k 

L   Aq   =  Aq 

n      * 
L   Aqn  =  Aq 

q     -  q   +  Aq 


■k  -k  -k 

where  Aq    and  Aq   are  intermediate  solutions.   The 
numerical  process  thus  results  in  a  block  tridiagonal  matrix 
inversion  for  each  coordinate.  The  size  of  the  blocks 
depends  on  the  number  of  elements  in  the  unknown  vector  q . 


Finite  Difference  Scheme  for  Axi symmetric 
Thin-Layer  Navier-Stokes  Equations 


In  a  manner  similar  to  the  three-dimensional  case,  the 
application  of  Beam  and  Warming  noniterative  implicit  scheme 
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to  the  axisymmetric  thin-layer  Navier-Stokes  equations  (19) 
results  in  the  following  finite  difference  equations: 


L,  L   Aq11  =  Qn  (27) 


'I    "C 


where 


s  At 

Qn  =    Aq11"1 (6rE  +  5  F+  6_G  -  Re_16_S)n 

1  +  s  1+s   *      *  g  g 

At  - 
Hn  -eF  J"1[(V?A,)2  +  (V_A_)2]Jqn 

1+s  *  *        ^  ^ 


and  L^  and  L£  are  given  in  equation  (26). 

Again,  the  Euler  implicit  scheme,  i.e.,  r=l  and  s=0,  is 
used  in  the  axisymmetric  thin-layer  Navier-Stokes  code  which 
results  in  first  order  accuracy  in  time  and  second  order 
accuracy  in  space. 


CHAPTER  IV 
THIN- LAYER  NAVIER- STOKES  CODES 


An  ax i symmetric  thin-layer  Navier-Stokes  code  and  a 
three-dimensional  thin-layer  Navier-Stokes  code  obtained 
from  the  U.S.  Army  Ballistic  Research  Laboratory  (BRL)  are 
employed  in  this  study  for  the  numerical  simulation  of 
transonic  turbulent  flows  past  a  projectile  model  with  a 
sting.  The  axi symmetric  code  is  a  product  of  BRL  and  NASA 
Ames  joint  effort  [11];  it  is  a  simplified  version  of  the 
three-dimensional  code  developed  at  NASA  Ames  [14].  The 
application  of  these  codes  to  transonic  projectile 
aerodynamic  problems  has  been  investigated  to  some  extent  by 
the  U.S.  Army  Ballistic  Research  Laboratory. 

The  three-dimensional  code  and  the  axi symmetric  code  are 
based  on  the  transformed  thin-layer  Navier-Stokes  equations 
(18)  and  (19),  respectively,  which  are  valid  only  for  high 
Reynolds  number  ideal  gas  flows.  In  these  codes,  the 
characteristic  velocity  used  in  the  transformation,  equation 
(8),  is  the  speed  of  sound  at  freestream  condition  c^ 
instead  of  the  freestream  velocity  V^;  consequently,  the 
freestream  Mach  number  M   in  equations  (10)  and  (12)  becomes 
unity.  The  solution  algorithm  adopted  for  the  transformed 
governing,  equations  (18)  and  (19)  is  a  noniterative, 
implicit,  approximate  factorization,  finite  difference 
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scheme  of  Beam  and  Warming  which  has  been  discussed  in  the 
previous  chapter.  It  is  noted  that  both  codes  have  an  option 
for  solving  inviscid  flow  problems  while  a  steady  state 
solution  is  obtained  as  the  converged  solution  of  the 
unsteady  impulsive  flow  problem. 

Both  the  axi symmetric  code  and  the  three-dimensional  code 
have  been  studied  in  detail.  Necessary  modifications  to  the 
computation  of  turbulent  eddy  viscosity  have  been  made  to 
the  axisymmetric  code  and  are  discussed  in  the  following; 
furthermore,  the  code  has  been  vectorized  for  Cray  X-MP 
supercomputers  to  improve  the  efficiency  of  computation. 
Moreover,  essential  subroutines  for  the  computation  of 
aerodynamic  forces  have  been  written  and  implemented  into 
both  axisymmetric  and  three-dimensional  codes. 

For  a  successful  application  and/or  modification  of  a 
Navier-Stokes  code,  one  must  be  familiar  with  certain 
relevant  features  of  the  code.  The  essential  features  of  the 
Navier-Stokes  codes  considered  and  the  subroutines 
implemented  are  described  and  discussed  in  the  following 
subsections  for  the  transonic  projectile  aerodynamics 
problems  investigated. 

Flow  Domain 

In  this  study,  flows  past  an  axisymmetric  projectile 
model  with  a  sting  are  considered  for  the  investigation.  For 
axisymmetric  flow  problems,  a  schematic  planar  flow  domain 


-  33  - 
and  the  corresponding  computational  domain,  ABCD,  are  shown 
in  figure  2.   Here,  the  line  AB  is  the  solid  boundary  of  the 
projectile  model,  the  line  BC  is  the  exit  or  downstream 
boundary  and  the  line  CD  is  the  outer  boundary  while  the 
line  AD  is  called  the  upstream  center  line.  For  a  projectile 
model  at  an  angle  of  attack,  that  is,  a  three-dimensional 
flow  problem,  the  flow  domain  in  the  physical  space  is 
generated  by  the  revolution  of  the  planar  domain  about  the 
axis  of  the  projectile  model.  Hence,  a  grid  network  for  a 
three-dimensional  problem  can  be  formed  by  a  sequence  of 
positions  of  the  planar  grid  ABCD  around  the  axis  of  the 
projectile  model.  The  corresponding  grid  system  in  the 
computational  domain  is  then  formed  by  stacking  the 
rectangular  grid  domain  ABCD.  Accordingly,  as  shown  in 
figure  2,  the  plane  ABB 'A'  is  the  solid  boundary  of  the 
projectile  model,  the  plane  BCC'B'  is  the  downstream 
boundary,  the  plane  CDD'C'  is  the  outer  boundary,  while  the 
plane  DAA'D1  is  a  singular  plane  representing  the  upstream 
center  line  AD. 

Boundary  Conditions 

For  a  selected  boundary-fitted  transformation,  such  as 
the  one  discussed  above,  suitable  boundary  conditions  must 
be  specified.   The  boundary  conditions  in  the  thin- layer 
Navier-Stokes  codes  are  applied  explicitly  and  require  the 
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specification  of  all  variables  p,  u,  v,  w  and  E  in  unknown 
vector  q  on  the  line  ABCD. 

The  application  of  solid  boundary  conditions,  line  AB,  is 
simplified  through  the  boundary  fitted  curvilinear 
coordinate  transformation  since  the  body  surface  has  been 
mapped  onto  the  plane  of  £=0.  The  no- slip  boundary  condition 
for  viscous  flow  is  enforced  by  setting  the  contravariant 
velocities  to  zero,  i.e.,  U=V=W=0  on  the  projectile  surface. 
For  inviscid  flow  problems,  a  linear  extrapolation  of 
tangency  contravariant  velocities  U  and  V  is  used  while  W=0 
for  impermeable  wall.   The  velocity  components  u,  v  and  w 
are  then  obtained  by  inverting  equations  (17).  Densities  at 
the  surface  are  assumed  to  be  equal  to  the  values  calculated 
on  the  first  grid  line  adjacent  to  the  boundary,  which  was 
always  well  within  the  viscous  sublayer.   The  pressure  along 
the  body  surface  is  obtained  from  the  three  transformed 
inviscid  momentum  equations.  By  combining  [ £  x ( £ -momentum)  + 

A 

£x  (  ti-  momentum)  +  £  x  (  ^-momentum)  ]  with  the  use  of  the 
continuity  equation,  metric  identities  (15)  and  neglecting 
viscous  terms,  one  finds 


11    L      +    £    C      +    t    C    )p,    +    (ti£      +ti£      +tiC)p 
v^xsx         y  y       ^zsz'^£        v  'x^x         y  y        'z^z'^n 

+     (r    2    +    r    2    +    r    2)p, 

X  Y  Z         <  (28) 

■     PfMt     +     UMx     +     V3lS     +     W3T<Z> 

"    Pu<«xu«    +    «yv*    +    CzV    -    PV^xut,    +    ^yvn    +    'zV 
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Since  viscous  terms  have  negligible  effect  on  surface 
pressure,  the  same  relation  (28)  has  been  used  in  viscous 
flow  computation  [14,16].   Equation  (28)  is  solved  at  the 
surface  by  using  second  order  central  differences  in  the  £ 
and  ti  directions  and  second  order  forward  differences  in  the 
C  direction.  The  total  energy  E  is  then  obtained  from  the 
known  pressure,  velocities  and  density  at  the  surface  by 
using  equation  (10). 

Along  the  upstream  line  of  symmetry,  line  AD,  which 
emanates  from  the  projectile  nose,  symmetry  conditions  are 
imposed  by  setting  the  second  order  finite  difference 
expression  for  the  first  derivative  equal  to  zero.  In 
three-dimensional  problems,  the  flow  variables  are 
determined  along  the  line  AD  independently  from  each  plane 
containing  line  AD.  The  averaged  values  obtained  from  all 
planes  are  then  used  on  the  singular  plane  DAA'D' .   Along 
the  far  field  outer  boundary,  line  DC,  freestream  conditions 
are  imposed.  On  the  downstream  boundary,  line  CB,  linear 
extrapolation  is  used  for  all  flow  variables  as  the  exit 
plane  flow  conditions  for  M  >1.   For  M  <1,  all  flow 
variables  are  obtained  from  linear  extrapolation  except  the 
total  energy  E  which  is  obtained  from  equation  (10)  by 
assuming  the  freestream  pressure  pw  fixed  at  the  exit  plane. 
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Initial  Condition 

A  steady  flow  problem  solved  by  the  Navier-Stokes  code  is 
treated  as  an  impulsive  flow  problem.   The  computational 
procedure  is  then  started  marching  in  time  until  an 
asymptotic  steady  state  is  reached.  Hence,  for  the 
projectile  flow  problem,  the  freestream  conditions  are  taken 
for  the  initial  conditions.   To  ease  the  sudden  jump  in 
solid  boundary  conditions  when  flow  impacts  the  projectile 
impulsively,  the  boundary  conditions  at  the  surface  are 
turned  on  slowly  over  the  first  30  time  steps.  Thus,  the 
boundary  conditions  of  the  unknown  vector  at  the  body 
surface  are  scaled  by 

q  =  qbc  x  SL  +  (1  -  SL)  q^ 


where 


sl  =  (io  -  15  x  r  +  6  x  r2)  r3 

r  =  (NC  -  l)/30 


where  NC  is  the  number  of  iterations,  q,  is  the  correct 
boundary  conditions  of  the  unknown  vector  which  has  been 
discussed  in  the  previous  section. 
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Viscosity 

The  effective  coefficient  of  viscosity  y  in  the 
Navier-Stokes  equations  is  comprised  of  two  parts,  the 
laminar  viscosity  yn  and  the  turbulent  eddy  viscosity  y^.. 
In  the  thin-layer  Navier-Stokes  codes,  the  evaluation  of 
laminar  viscosity  y»  as  a  function  of  temperature  is 
obtained  from  Sutherland's  theory  of  viscosity  and  given  by 


*  _  »«  .  .  T   .1.5  T-+  S  -  _.1.5  ,  !  +S"  , 

V       T        T  '+  S  T  +S 

oe  oo 


where  the  constant  S  is  198. 6°R  [15]  and  S  =3/1^. 

The  turbulent  eddy  viscosity  model  employed  in  the  codes 
is  the  one  proposed  by  Baldwin  and  Lomax  [2].  It  is  a 
two-layer  algebraic  model  in  which  an  eddy  viscosity  is 
calculated  for  an  inner  and  an  outer  region.  The  inner 
region  follows  the  Prandtl-Van  Driest  formulation 


dinner  "  P*2  I  •  I  (29> 


with    I   =  ky[l  -  exp(-y+/A+) 


where  y  is  the  normal  distance  to  the  surface  and  |w|  is  the 
magnitude  of  the  vorticity 


|w|  =  [(3yu-3xv)2  +  (3zv-3yw)2  +  (3xw-3zu)2]1/2 
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and  y   is  the  law  of  the  wall  coordinate 


~*     (PwTw>1/2y 

y  = 

"w 


where  t  ,  p   and  u   are  the  local  shear  stress,  density  and 
laminar  viscosity  evaluated  at  the  wall,  respectively. 

The  eddy  viscosity  for  the  outer  region  is  given  by  the 
Clauser  formulation 

Pouter  =  PKCcpFwakeFkleb(y>  <30> 


where  K  is  the  Clauser  constant,  C is  an  another  constant 

'   cp 

and  F   >e  is  the  smaller  value  of  the  following: 


wake  '  ^max  max 

wake  '   wk^maxudif  '  max 


The  quantities  y__„  and  F__„  are  determined  from  the 
function 


F(y)  =  y|w|[l  -  exp(-y+/A+)] 


where  F__v  is  the  maximum  value  of  F(y)  and  "y__„  is  the 

ITlclX  UlaX 

value  of  y  at  which  it  occurs.  The  Klebanoff  intermittency 
factor  Fkleb(y)  is  given  by 
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Fkleb(y)  ■  [1  +  S-^Ckleby/ymax)^"1 


and  the  quantity  ud^£  is  the  difference  between  maximum  and 
minimum  total  speed 

udif  =  (u*  ♦  v2  t  w^)maxV2.(u2  ,  v2  +  w2)minl/2 

where  the  minimum  total  speed  is  zero,  unless  the  spinning 
projectile  is  considered. 

Rather  than  attempting  to  match  (v^)inner  and  (v^)outer 
values  in  the  overlap  region,  both  (v^)inner  and  (v^)outer 
are  calculated  at  every  grid  point  in  the  flow  field  and  the 
smaller  value  is  used  as  turbulent  eddy  viscosity.  The 
constants  used  for  this  model  are  A   =  26,  C    =  1.6,  C  k  = 
0.25  and  Ckleb  =0.3  .   Note  that  in  this  section,  all 
quantities  are  dimensional  except  those  with  a  superscript 
* . 

In  the  original  axisymmetric  thin-layer  Navier-Stokes 
code,  the  turbulent  eddy  viscosity  was  calculated  only  for 
the  first  25  grid  points  off  the  projectile  surface  in  order 
to  increase  computational  efficiency.  However,  in  the 
current  grid  networks  used,  the  25th  point  is  well  within 
the  boundary  layer  and  subsequent  zero  values  for  the 
turbulent  eddy  viscosity  in  the  remaining  region  of  the 
boundary  layer  led  to  instabilities  in  the  numerical 
algorithm.  The  code  has  been  modified  to  include  the 
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~ turbulent  eddy  viscosity  at  every  grid  point  in  the  entire 
domain. 


Forces  and  Moments 

The  motion  of  the  fluid  relative  to  a  solid  boundary 
exerts  a  dynamic  force  on  the  surface  which  consists  of  both 
frictional  and  pressure  forces.  These  surface  forces  can 
create  the  pitching,  yawing  and  rolling  moments  relative  to 
the  center  of  gravity.   In  the  cases  we  considered,  the 
axi symmetric  projectile  without  spinning,  there  is  no 
contribution  to  the  side  (lateral)  force,  nor  the  yawing  and 
rolling  moments.  Note  that  both  pressure  and  shear  stresses 
can  cause  the  pitching  moment  for  the  projectile  at  an  angle 
of  attack. 

The  forces  and  the  pitching  moment  acting  on  an  arbitrary 
shape  can  be  obtained  by  integration  of  the  surface  pressure 
and  shear  stress  distributions  and  moments  about  the  center 
of  gravity,  respectively,  as  follows: 

(1)  Force  components  in  axial  direction 


FpA  =  SS   p  sin9  dA 
FfA  =   II    x    cos9  dA 


(31) 
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(2)  Force  components  in  normal  direction 


F  N  =   -    II   p   cos8  cos0  dA 
F^j,  =    II    x    sin8  cos0  dA 


(32) 


(3)  Pitching  moment  about  the  center  of  gravity 

MC.G.  =  U    [FN  (XC.G.  "  X)  +  FA  Z  cos0]  dA       (33) 

Here,  p  and  t  are  the  pressure  and  shear  stress 
distributions  on  the  body.  The  subscripts  p  and  f  refer  to 
the  forces  due  to  pressure  and  shear  stresses,  respectively, 
and  subscripts  A  and  N  indicate  the  force  components  in  the 
axial  and  normal  direction,  respectively.   The  angle  <j> 
starts  on  the  leeward  side  and  measures  in  the 
circumferential  direction  and  9  is  the  inclination  of  local 
slope  of  the  surface,  as  shown  in  figure  3 .  Xc  G   is  the 
axial  position  of  the  center  of  gravity  measured  from  the 
nose  of  projectile  which  is  3.6  for  this  particular 
projectile.   X  is  the  local  axial  position  and  Z  is  local 
position  in  the  radial  direction,  measures  from  the  axis  of 
the  projectile.   The  pitching  moment  is  considered  positive 
when  it  is  in  the  nose-up  sense. 

With  the  same  dimensionless  quantities  in  equations  (8) 
and  characteristic  velocity  cm,    forces  in  the  equations  (31) 
and  (32)  can  be  put  into 
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?      ?   *       * 
F   =Lp   c    F    =6F 

H  P       F  (34) 

Ff  =  L  »».  c~  Ff*  =  P  —  Ff* 

Re 


The  ratio  of  the  forces  due  to  pressure  and  shear  stress  is 
therefore  F    to  Ff  /Re. 

The  pressure  and  shear  stress  coefficients  are  expressed 
as 


*      * 
P  -  P„     P  "  P. 

C 


d/2)p  V  2    (1/2)M  2 


* 

t  /Re 


(35) 


C*  = 


(1/2)P  V  2    (1/2)M  2 


and  the  pitching  moment  coefficient  as 


MC.G.         MC.G.* 
c   _  ^^ _  ^^ ^36^ 

(V2)PV2L    (1/2)M  2L 


Because  the  moment  has  dimensions  of  a  force-length  product, 
the  additional  length  scale  L  in  the  denominator  of  moment 
coefficient  is  introduced. 

The  component  of  the  resultant  force  in  the  direction  of 
the  freestream  velocity  V   that  passes  the  body  is  the  drag, 
which  is  one  of  the  most  important  parameters  in  aerodynamic 
designs.  It  can  be  expressed  as 
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FD  =  FN  sina  +  FA  COSa  (37) 

where  a  is  the  angle  of  attack.  The  corresponding  lift  force 
is 

FL  =  FN  cosa  -  FA  sina  (38) 

For  axi symmetric  flow  cases,  i.e.  a=0°,  the  drag  force  will 
be  the  axial  force  and  lift  force  will  be  zero. 

The  computation  of  forces,  moment  and  shear  stress 
distribution  has  been  added  to  both  thin-layer  Navier-Stokes 
codes.  Thus,  the  computer  codes  have  become  more  complete  in 
evaluating  the  flow  features. 


Rate  of  Convergence 

The  numerical  algorithm  is  a  time-marching  scheme  in 
which  a  steady  state  solution  is  obtained  as  a  time 
asymptotic  solution  of  the  unsteady  equations.  The 
convergence  of  the  process  can  be  observed  from  the  average 
residual . 


(T    B      2.1/2 

Average  residual  =  J  —   (39) 

A  t  •  (  JMAX-  2  )  ( KMAX- 2  )  ( LMAX-  2  ) 
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where  R  is  either  the  vector  of  RHS(26)  or  RHS(27).   The 
total  numbers  of  grid  points  in  the  £,  n  and  £  direction  are 
JMAX,  KMAX  and  LMAX,  respectively.  And  summations  over  /,  k, 
I   range  from  2  to  JMAX-1,  KMAX-1,  LMAX-1,  respectively.  The 
summation  of  m  is  from  1  to  the  total  number  elements  of 
unknown  vector  q,  which  for  this  case  is  5. 

For  steady  state  solutions,  from  equation  (26)  or  (27),  we 
know  that  the  vector  R  is  zero  since  the  solution  will  no 
longer  vary.  This  implies  then  that  the  average  residual 
becomes  zero  when  the  steady  state  solution  is  reached.  In 
the  numerical  process  of  this  study,  the  convergence 

criterion  for  the  steady  state  solutions  being  reached  has 

-4 
been  set  to  the  average  residual  less  than  10 


CHAPTER  V 
HYPERBOLIC  GRID  GENERATION 


In  the  applications  of  thin-layer  Navier-Stokes  codes,  a 
grid  network  must  be  provided.  The  grid  network  used  in  this 
study  is  generated  from  a  grid  generation  code,  named 
GRIDGEN,  obtained  from  the  U.S.  Army  Ballistic  Research 
Laboratory.   This  code  can  provide  a  grid  network  generated 
either  by  elliptic  or  by  hyperbolic  equations.   For  external 
flow  problems,  hyperbolic  grid  generation  is  feasible. 
Since  specifying  the  precise  positions  and  shape  of  the 
outer  boundary  is  not  important,  the  grid  around  the 
geometry  of  interest  can  be  obtained  by  the  specification  of 
grid  positions  at  one  boundary,  and  then  marched  outward 
from  this  boundary  until  the  far  field  is  reached.   The 
hyperbolic  grid  generation  is  a  noniterative  scheme,  so  the 
required  computer  time  is  almost  equivalent  to  that  of  one 
iterative  sweep  in  solving  elliptic  grid  generation 
equations.  Furthermore,  this  hyperbolic  grid  has  the 
orthogonal  property,  and  grid  smoothness  can  be  controlled 
by  a  selected  grid  distribution  function. 
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Governing  Equations 

Construction  of  hyperbolic  planar  grid  (x,z)  proceeds  by 
selecting  two  differential  equations  which  serve  as 
constraint  equations  that  are  solved  to  obtain  the 
coordinates  of  grid  points  [17].   It  is  known  [8]  that  an 
orthogonal  grid  will  reduce  the  truncation  error  terms 
associated  with  finite  difference  approximation  formulated 
in  generalized  coordinates;  thus  one  equation  selected  is 
that  which  enforces  orthogonality 


U  •  V£  =  0      or      £XCX  +  KZCZ   =   0  (40) 


For  a  two-dimensional  general  coordinate  transformation, 
the  metric  relations  and  Jacobian  J  of  the  transformation 
are 


z_                                                           X^  Z£  X>- 

Kx  =  —        Kz  =  -  —  cx  =  -  —  cz  =  —                (41) 

XJ         ZJ  XJ  ZJ 

J   =   x^    -    x^  (42) 


The  equation  (40)  can  then  be  written  in  the  computational 
space  ( C / C )  as 

X;X   +  ZfZ_  =  0  (43) 
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The  second  equation  is  chosen  to  specify  a  finite  grid 
cell  area  to  insure  a  nonsingular  mapping  through  the 
transformation.   In  numerical  implementations,  A5=A£=1,  the 
Jacobian  in  equation  (42)  is  approximately  equal  to  the 
physical  cell  area.   A  set  of  hyperbolic  grid  generation 
equations  is  then  formulated  by  equations  (43)  and  (42). 

Equations  (43)  and  (42)  are  nonlinear  equations  which  can 
be  locally  linearized  by  expanding  about  a  known  state 
(x0,z°).   Neglecting  the  high  order  term,  e.g., 
(x-x° ) £ (x-x° )  ,  we  have 

x^x^  =  [x°  +  (x  -  x0)]^  [x°  +  (x  -  x°)]? 
=  x?°X£  +  x^°x?  -  x^xc° 

etc.  where  the  superscript  °  denote  the  known  state.  Upon 
substituting  these  relations  into  equations  (43)  and  (42), 
one  finds 


x  °xf  +  z  °z>.  +  x,°x   +  zf"z   =  0 

6   *     C   *     *  Q  K       C  (44) 

z  °X£  -  x  °Z£  -  z/x   +  x>'z   =  J  +  J° 


where  the  grid  areas  J  are  to  be  specified  and  J°  is 
determined  at  the  known  state.  Therefore,  all  quantities  of 
RHS(44)  are  known.  Equations  (44)  can  be  put  into  the 
compact  form 
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AW,    +   BW     =  f 


(45) 


with 


A  = 


x  °   z  ° 

Z   °   -X   ° 


B  = 


X   °   7   ° 
-■7   °   V   ■ 


f  =\  I  w  = 

|j  +  J° 


Assume  that  the  £  coordinate  runs  in  the  direction  away 
from  a  specified  boundary.  Then  the  governing  equations  for 
the  propagation  problem  can  be  written  as 


where 


CW^    +  Wc  =  g 


C   =  B~1A 


g  =  B_1f 


The  eigenvalues  of  C  are  two  distinct  real  roots  which 
implies  that  the  governing  equations  (45)  form  a  hyperbolic 
system. 


Finite  Difference  Form 


In  the  GRIDGEN  code,  the  finite  difference  approximations 
used  for  the  hyperbolic  grid  equations  (45)  are  second  order 
central  difference  for  £  derivatives  and  first  order 
backward  difference  approximation  for  derivatives  in  the 
marching  direction  C;  i.e., 
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Then,  an  implicit  finite  difference  scheme  can  be  formed 


1  1 


2 

(46) 


which  forms  a  2x2  block-tridiagonal  set  of  equations  to  be 
solved  on  each  successive  constant  c-line.   Note  that  the 
column  matrix  fy  £+j  on  RHS(46)  is  assumed  known,  since  it 
only  involves  the  quantities  J+J° .   The  grid  areas  J  on  the 
next  Z+1    constant  {-line  can  be  specified  in  any  way.   In 
the  GRIDGEN  code,  it  has  been  specified  as 


1 
Jj,l+1    =  ~   AL  x  ASi 


where  AL  is  the  total  length  of  two  adjacent  cells,  as  shown 
in  figure  4,  and  AS£  is  the  specified  height  of  two 
successive  constant  C-lines.  The  distribution  AS  used  in  the 
code  is  an  exponential  distribution,  and  is  determined  by 
the  relation 
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Figure  4.  Spec 


ification  of  grid  area  in  GRIDGEN  code. 
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AS-t  =  AS0  (1  +  e^~1  '    *•   =    1/IMAX-l       (47) 


where  IMAX  is  the  total  number  of  grid  points  along  the 
given  arc  length,  and  ASQ  is  the  specified  grid  spacing  at 
one  end  of  arc.  The  parameter  e  is  determined  by  a 
Newton-Raphson  iteration  process  so  that  the  sum  of  the 
above  increments  matches  the  known  arc  length.   The  grid 
spacing  ASQ  next  to  the  boundary  considered  in  this  study  is 
0.01  for  inviscid  flow  problems;  however,  for  viscous 
turbulence  cases,  first  grid  spacing  ASQ  of  2xl0-5  is 
required. 

Several  different  distribution  functions  have  been 
discussed  and  investigated  in  detail  [19,20,22],  and,  in 
comparison,  a  hyperbolic  tangent  distribution  has  the  better 
overall  distribution  [19,20].   For  extensive  study  of  the 
effectiveness  of  the  grid  network  on  the  flow  solution,  a 
hyperbolic  tangent  distribution  is  implemented  in  the 
computer  code.   A  hyperbolic  tangent  distribution  is 
constructed  by 


6    -t  -  1 

tanh[-( 1) 

2  IMAX  -  1 

Si   =  Stot  *1  +  ) 

tanh(5/2) 

(48) 

LSi   =   Si+1    ~    Sl  ,    i   =    1.     IMAX-1 
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where  Stot  is  the  total  length  of  the  arc  and  the  parameter 
6  is  determined  by 


=  (IMAX  -  1) 


ASQ 


sinh(5)  stQt 

Two  computational  planar  grids  for  the  projectile 
configuration  with  the  same  boundary  distribution  and 
ASq=2x10"''  grid  but  different  normal  distributions, 
equations  (47)  and  (48),  for  the  projectile  configuration 
are  shown  in  figure  5 .   We  observed  that  the  grid  with 
hyperbolic  tangent  clustering  provides  a  smoother  grid  away 
from  the  projectile  than  the  one  with  exponential  clustering 
in  the  radial  C-direction.  However,  near  the  body  surface, 
the  grid  with  exponential  clustering  is  smoother  than 
hyperbolic  clustering. 

A  hyperbolic  grid  generation  is  a  propagation  problem;  it 
requires  the  initial  grid  points.  The  grid  point 
distribution  on  the  body  surface,  in  this  study,  is  adapted 
to  the  flow  solution  in  the  £ -direction  and  then  propagated 
outward  in  ^-direction.  Furthermore,  a  three  dimensional 
grid  for  the  axi symmetric  projectile  geometry  considered  is 
generated  by  rotating  a  two  dimensional  planar  grid  about 
the  axis  of  symmetry  and  maintaining  a  uniform  angular  space 
between  the  rotated  planes. 
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(a)  Based  on  exponential  clustering. 
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(b)  Based  on  hyperbolic  tangent  clustering. 

Figure  5.  90x40  hyperbolic  grids  with  different  clusterings 
in  the  normal  direction. 


CHAPTER  VI 
FLOW  PROBLEM 


The  objectives  of  this  study  are  to  gain  further 
understanding  of  complex  flow  simulations  using  the 
thin- layer  Navier-Stokes  approximation  and  to  investigate 
the  relative  importance  of  the  viscous  force  to  the  pressure 
force  on  a  projectile  at  transonic  speed. 

The  geometric  configuration  considered  is  an  axisymmetric 
secant-ogive-cylinder-boattail  (SOCBT)  projectile.   The 
projectile  model  has  a  3-caliber  secant-ogive  part  followed 
by  a  2-caliber  cylinder  and  1-caliber  7-degree  boattail  as 
shown  in  figure  6.  The  model  was  supported  by  a  base-mounted 
sting  in  the  experiments.  Experimental  data  are  available 
for  the  pressure  distribution  along  the  projectile  surface 
for  a  range  of  Mach  numbers  in  the  transonic  regime,  0.91, 
0.94,  0.96,  0.98,  1.10  and  1.20;  at  angles  of  attack  of  0°, 
2,4,6   and  10°;  and  at  circumferential  positions  around 
the  model  in  22.5°  increments  [7].  The  computational  model 
of  the  projectile,  however,  is  formed  by  extending  the 
boattail  part  of  the  projectile  for  another  1.77  calibers  to 
meet  the  sting.  This  modification  is  made  to  avoid  the 
difficulty  of  simulating  the  base  flow  region. 

For  the  external  flow  problems,  which  are  considered  here, 
the  outer  boundary  needs  to  be  specified.  Theoretically,  the 
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outer  boundary  should  be  as  far  as  possible  from  the 
projectile,  so  that  freestream  conditions  can  be  imposed. 
The  wind  tunnel  experiments  supply  the  freestream  stagnation 
temperature  T0=580°R;  moreover,  the  projectile  body  is 
assumed  adiabatic  and  impermeable. 


CHAPTER  VII 
RESULTS  AND  DISCUSSION:  AX I SYMMETRIC  FLOWS 


For  the  purposes  of  evaluating  the  accuracy  and 
applicability  of  the  thin-layer  Navier-Stokes  codes 
described  previously  for  complex  flow  problems,  comparisons 
of  measured  and  predicted  aerodynamic  characteristics  and 
associated  flow  fields  for  various  flow  conditions  have  been 
made.  The  finite  difference  scheme  for  the  Navier-Stokes 
equations  requires  the  additional  artificial  second  order 
implicit  and  fourth  order  explicit  dissipation  terms, 
equations  (24)  and  (25)  respectively,  to  control  numerical 
instabilities.  According  to  a  linear  stability  analysis,  the 
scheme  is  unconditionally  stable  for  linear  equations  when 
eI=2eE  I3'4'2*]-   This  weight  is  used  in  all  subsequent 
investigation  of  axisymmetric  flow  problems. 


Inviscid  Flow 
Flow  Domain  Determination 

For  SOCBT  projectile  model,  the  experimental  data  are 
available  in  a  range  of  M^O.91  to  1.20.  A  flow  case  of 
M«,=0-96  is  selected  for  detailed  investigation.   In  order  to 
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determine  adequate  flow  domains,  the  axisymmetric 
Navier-Stokes  equations  (19)  were  solved  on  a  90x40 
hyperbolic  grid  with  exponential  clustering,  equation  (47), 
in  the  normal  direction.  The  two  domains  used,  St  t=18  and 
24  are  about  3  and  4  times  the  projectile-length, 
respectively.  The  results  presented  in  figure  7  show  that 
the  computed  surface  pressures  are  very  similar  and  give  an 
accurate  prediction  of  the  experimental  data  except  on  the 
boattail  region.   The  results  indicate  that  both  domains  18 
and  24  seem  to  be  sufficiently  large  for  imposing  freestream 
conditions  on  the  outer  boundaries.  Nevertheless,  in  theory, 
the  better  solution  is  obtained  with  a  larger  domain  when 
enough  grid  resolution  is  provided. 

Effect  of  Artificial  Dissipations 

Since  the  flow  domain  of  18  can  provide  accurate  results, 
the  78x28  hyperbolic  grid  with  domain  18  was  used  with  two 
different  sets  of  dissipation  terms:  one  is  eI=2eE=4At; 
another  is  eI=2EE=8At  for  the  inviscid  flow  case.  The 
influence  of  artificial  dissipation  terms  on  the  flow 
solutions  is  presented  in  figure  8.  When  more  dissipation  is 
added,  smearing  of  the  numerical  solution  occurs  near  the 
shock/boundary  layer  interaction'  region,  about  at  the  middle 
of  cylinder  part,  since  the  dissipative  nature  tends  to 
smear  the  predicted  shock  over  several  mesh  intervals.  In 
order  to  prevent  the  deterioration  of  the  solution,  the 


M  =  .96 


i 

with  a  90x40  grid  but  different  domain. 


Stot=18 


Stot=24 


M  -  .36 


Figure  8.  Computed  surface  pressure  of  inviscid  solution 
with  a  78x28  grid  but  different  dissipations 
:  Ej=2EE=4At  


Ej=2eE=8At 
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added  dissipation  terms  should  be  as  small  as  possible; 
however,  the  solution  algorithm  became  unstable  when 
Ej=2eg=2At  is  employed. 


Viscous  Flow 

In  turbulent  flow  applications,  the  closure  turbulence 
model  used  in  the  codes  is  a  two- layer  algebraic  eddy 
viscosity  turbulence  model  of  Baldwin  and  Lomax.   The 
computed  results  obtained  show  that  a  90x60  hyperbolic  grid 
with  exponential  clustering  in  ^-direction  and  flow  domain 
Stot=24,  which  is  about  4  times  the  projectile  length,  can 
give  very  accurate  solutions.  Figure  9a  shows  that  the 
computed  surface  pressure  coefficient  C   is  in  excellent 
agreement  with  measured  data  while  the  computed  shear  stress 
T££  presented  in  figure  9b  indicates  that  there  is  no  flow 
separation  on  the  projectile  surface.  The  corresponding 
pressure  and  Mach  contours  are  shown  in  figure  10.  The  dense 
regions  in  the  contour  plots  indicate  that  large  gradients 
occur;  i.e.,  it  implies  the  corresponding  flow  property 
changes  rapidly.  The  coalescence  of  the  Mach  lines  to  the 
right  of  the  supersonic  region  represents  the  position  of 
the  shock.  Two  shock  locations,  one  in  the  middle  of 
cylinder  part  and  another  around  the  middle  of  boattail 
portion,  agree  quite  well  with  a  shadowgraph  in  reference  7. 
The  expansion  waves  also  are  shown  at  the  ogive-cylinder  and 


-  63  - 


-.4- 


M  =    .96 


(a).  Pressure  coefficient: 


0:  measured  data. 


(b).  Surface  shear  stress  t 


U 


Figure  9.  Computed  results  with  a  90x60  grid. 
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(a)    Pressure   contours 
p:    0.48->0.86,    Ap=0.02 


8  2  4  e 

(b)    Mach   contours 
M:    0.88^1.0,     AM=0.01;    M:     1.0-»1.26,     AM=0.02 


Figure    10.    Contour  plots   of  1^=0.96. 
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cylinder-boattail  junctures.  The  computed  nondimensional 
drag  force,  D=FDxlO~  ,  due  to  surface  pressure  is  D  =40.28 
over  the  6-caliber  projectile,  excluding  the  base  drag, 
while  the  drag  force  resulting  from  skin-friction  is 
Df=20.53.  This  shows  the  importance  of  viscous  flow 
computation  for  accurate  aerodynamic  force  prediction. 

Effect  of  Artificial  Dissipations 

The  effect  of  artificial  dissipation  terms  upon  the 
accuracy  of  numerical  solutions  of  turbulent  flow  cases  was 
next  considered.  The  same  grid  network  90x60  of  domain 
Stot=24  was  used  but  with  eI=2eE=8At.  The  computed  C   as 
shown  in  figure  11a  seems  to  agree  rather  well  with  the 
accurate  results,  while  the  quality  of  the  solution  is 
degraded  near  the  shock/boundary  layer  interaction  region 
owing  to  the  nature  of  dissipation.   The  resulting  pressure 
drag  is  Dp=32.91  which  is  18%  less  than  that  of  the  accurate 
solution.  Figure  lib  shows  the  difference  between  the 
calculated  shear  stress  distributions;  however,  the 
resulting  shear  drag  is  about  the  same,  Df=20.12.   The 
difference  in  the  drag  forces  due  to  pressure  is  mainly  due 
to  the  difference  in  the  predicted  pressure  on  the  boattail 
part,  since  the  pressure  intensity  on  the  cylinder  portion 
does  not  contribute  to  the  drag  force  in  axi symmetric  flow 
cases. 
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(a).  Pressure  coefficient;     0:  measured  data 


(b).  Surface  shear  stress  !_► 

Figure  11.  Computed  results  with  a  90x60  grid  but 
different  dissipations. 


ej=2EE=4At 


z j=2eE=8At 
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Effect  of  Eddy  Viscosity  Distribution 

The  effect  of  an  averaged- technique  originally 
implemented  in  the  axi symmetric  thin- layer  Navier-Stokes 
code  for  computing  eddy  viscosity  also  has  been  investigated 
on  90x60  grid  at  M<x>=0.96.  Figure  12  shows  the  distribution 
of  eddy  viscosity  with  and  without  the  averaging  technique 
at  three  different  boundary  point  stations  identified  in 
figure  11a.  It  is  clear  that  the  averaging  scheme  yields 
improper  zig-zag  eddy  viscosity  distribution  above  the 
viscous  sublayer;  nevertheless,  it  has  negligible  effect  on 
surface  pressure.   In  fact  the  computed  C  -distribution  is 
almost  exactly  the  same  as  that  of  the  accurate  solution 
shown  in  figure  13;  however,  the  shear  stress  computed  is 
consistently  less  than  that  of  the  accurate  solution  over 
the  entire  projectile  surface.  The  corresponding  pressure 
drag  force  is  D  =41 . 41  and  skin-friction  drag  is  Df=18.99, 
which  are  different  by  3%  and  8%  respectively. 

Effect  of  Grid  Resolution  in  Streamwise  Direction 

The  grid  networks  78x40  and  90x40  with  domain  St  t=18  are 
used  for  comparison.  Of  the  12  additional  points  in  the 
^-direction  of  the  90x40  grid  network,  2  points  were  added 
in  the  ogive  part  and  10  points  were  added  in  cylinder  part. 
From  the  results  presented  in  figure  14,  we  observed  that 
the  solutions  are  generally  quite  the  same  except  in  the 
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(a)  Sample  eddy  viscosity  distributions  based  on 
an  averaged- scheme  implemented  in  the  original 
axisymmetric  thin-layer  Navier-Stokes  code 
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Figure  12.  (b)  Sample  eddy  viscosity  distributions  computed 
without  averaging. 
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(a).  Pressure  coefficient;      Q:  measured  data. 


(b) .  Surface  shear  stress  t 
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Figure  13.  Computed  results  with  and  without  an 
averaged- scheme . 
.:  without  averaged- scheme    .  with  averaged- scheme 
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6.  XSO 


Figure  14.  Computed  surface  pressure  with  domain  S+  <_=18 
but  different  grid  resolution  in  £-direc?ion. 
•"  78x40  grid     ■  90x40  grid. 
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cylinder  portion  which  is  a  shock/boundary  layer  interaction 
region.  In  shock/boundary  layer  interaction  regions,  the 
solutions  are  more  sensitive  to  grid  resolution.  The  results 
show  that  the  90x40  grid,  with  more  points  in  the  cylinder 
portion,  has  indeed  improved  the  solution  characteristics, 
especially  right  after  the  ogive-cylinder  juncture.  However, 
the  computed  C  -distribution  of  both  cases  do  not  reach  the 
peak  value  of  experimental  data  on  the  cylinder  part.  This 
could  be  attributed  to  insufficient  domain  or  grid 
resolution  on  the  shock/boundary  layer  interaction  region, 
etc. 

Effect  of  Grid  Resolution  in  Normal  Direction 

To  further  test  for  the  required  grid  resolution  in  the 
normal  direction,  the  90x40  grid  with  domain  StQt=24  was 
solved.  The  results  obtained  and  shown  in  figure  15  are 
compared  with  the  most  accurate  solution  from  the  90x60 
grid.   Both  grid  networks  have  the  same  grid  distribution  on 
the  body  surface  and  the  same  exponential  clustering 
function  but  different  numbers  of  grid  points,  40  and  60,  in 
the  C-direction.  As  is  evident  in  figure  15  by  the  poor 
accuracy  on  the  cylinder  and  boattail  portions,  the  90x40 
grid  did  not  provide  enough  resolution.   The  computed  drag 
force  due  to  pressure  is  23.95  which  is  about  41%  less  than 
that  of  the  90x60  case.  However,  the  drag  force  due  to  skin 
friction  is  20.28  which  has  no  significant  difference  from 
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(a).  Pressure  coefficient;     0:  measured  data 
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(b) .  Surface  shear  stress  x 


Figure  15.  Computed  results  with  St  .=24  but  different 
grid  resolution  in  C-direction. 


:  90x40  grid 


U 
bu 

90x60  grid. 


-  73  - 
20.53  of  the  90x60  case,  even  though  the  shear  stress 
distributions  are  quite  different,  as  shown  in  figure  15b. 

The  flow  problem  with  domain  Stot=24  has  been  solved  again 
with  a  90x40  hyperbolic  grid  based  on  hyperbolic  tangent 
clustering,  equations  (48),  in  the  normal  direction.   The 
characteristics  of  exponential  and  hyperbolic  tangent 
clustering  functions  have  been  discussed  in  chapter  V;  the 
corresponding  grid  networks  are  shown  in  figure  5.   The 
computed  Cp  presented  in  figure  16a  shows  that  the  grid 
resolution  using  40  points  with  exponential  clustering  in 
the  normal  direction  is  not  sufficient;  moreover,  a  grid 
which  resulted  from  the  use  of  hyperbolic  tangent  clustering 
yields  a  better  result;  it  yielded  D  =27.38  while  the  other 
grid  resulted  in  D  =23.95.  These  pressure  drag  forces  are, 
respectively,  32%  and  41%  less  than  the  accurate  value  of 
Dp=40.28.  The  computed  shear  stress  given  in  figure  16b 
indicates  that  the  shear  stress  distribution  is  very 
sensitive  to  the  pressure  field,  particularly  in  the 
shock/boundary  layer  interaction  region;  however,  the 
resulting  shear  drag  forces  are  about  the  same,  20.05  and 
20.28,  respectively.  It  should  be  pointed  out  that  a  grid 
with  hyperbolic  tangent  clustering  does  not  always  give  more 
accurate  results.  In  fact  the  flow  problem  with  a  domain 
Stot=18  also  has  been  solved  on  the  78x28  grids.  The  results 
obtained  and  presented  in  figure  17  show  that  a  grid  with 
hyperbolic  tangent  clustering  gives  better  result  on  the 
boattail  portion  but  the  other  grid  with  exponential 
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(a).  Pressure  coefficient; 


0  :  measured  data. 


(b) .  Surface  shear  stress  t 
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Figure  16.  Computed  results  with  a  90x40  grid  but  different 
clustering  in  C-direction. 

:  with  exponential  clustering 

:  with  hyperbolic  tangent  clustering. 
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M  =  .96 
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Figure  17.  Computed  results  with  a  78x28  grid  but  different 
clustering  in  ^-direction. 

-- — -  :  with  hyperbolic  tangent  clustering. 
— •  :  with  exponential  clustering 
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clustering  yields  a  much  more  accurate  solution  on  the 
cylinder  portion  of  the  projectile. 

Importance  of  Viscous  Drag  to  Pressure  Drag 

In  order  to  gain  further  insight  as  to  the  relative 
importance  of  the  viscous  drag  to  the  pressure  drag,  the 
sequence  of  transonic  flows  1^=0.91,  0.94,  0.98,  1.10  and 
1.20  past  the  SOCBT  projectile  model  at  zero  angle  of  attack 
have  also  been  solved  with  the  same  90x60  hyperbolic  grid 
using  exponential  clustering  in  the  normal  direction.  We 
observed  by  comparing  computed  C  -distribution  with 
experimental  data,  shown  in  figures  18,  19,  20,  21  and  22, 
that  the  grid  which  is  a  very  good  grid  for  M  =0.96 
apparently  is  not  the  one  best  suited  for  the  other  flow 
cases;  however,  the  results  computed  are  acceptably 
accurate.  The  shear  stress  distributions  computed  have  shown 
that  there  is  a  small  region  of  the  reversed  flow  at  the 
boattail  portion  for  M^O.91  and  0.94;  surprisingly  the 
computed  C   on  the  boattail  portion  is  in  excellent 
agreement  with  measured  data.  It  seems  to  imply  that  the 
Baldwin-Lomax  turbulence  model  can  give  an  accurate  solution 
for  small  separated  flow  regions.  The  drag  forces  due  to 
pressure  and  shear  stress  are  listed  in  the  following  table. 
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Figure  18.  Computed  results  with  a  90x60  grid. 
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M  =    .94 


(a).  Pressure  coefficient; 


0:  measured  data, 
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Figure  19.  Computed  results  with  a  90x60  grid. 
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M  =    .98 


'"""I" I'" |MM..I.I, ■■■..,,■.,,      ,,, , 


2  4 

(a).    Pressure   coeffici 


6       X'D 

ent;  0   :    measured  data, 


M  =    .98 


■ 1 1 1 1 1 1| i  ii  M 1 1 1 1 1 1 1 1 1 1 1 1 1  n 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1  r 1 1 1 1  »i  i 


■  1 1 1 1 1 1 1 1 1 


2  4  6   X/D 

(b).  Surface  shear  stress  t 


Figure  20.  computed  results  with  a  90x60  grid. 
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Figure  21.  Computed  result 


s  with  a  90x60  grid. 
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M  =    1.2 


(a).  Pressure  coefficient; 
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(b).  Surface  shear  stress  t 
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Figure  22.  Computed  results  with 


a  90x60  grid. 
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Table  2.  Drag  forces  due  to  pressure  and  shear  stress 
at  different  Mach  numbers. 


MM     0.91     0.94     0.96     0.98     1.10     1.20 


Dp    20.81    32.22    40.28    63.78    78.95    98.24 
Df    18.35    19.33    20.53    21.88    27.67    32.56 


The  data  in  the  table  indicate  that  the "forces  due  to 
pressure  as  well  as  due  to  skin  friction  are  increased  with 
increasing  Mach  number.  The  ratio  of  the  computed  pressure 
drag  Dp  to  skin  friction  drag  Df  is,  respectively,  1.13, 
1.67,  1.96,  2.92,  2.85  and  3.02  for  flow  cases  M  =0.91, 

oo  ' 

0.94,  0.96,  0.98,  1.10  and  1.20.  These  computed  data 
indicate  that  the  skin  friction  drag  is  rather  important  in 
the  low  transonic  regime,  and  is  not  as  significant  in  the 
higher  transonic  regime. 

The  corresponding  pressure  and  Mach  contour  plots  of 
M^O.91,  0.94,  0.98,  1.10  and  1.20  are  shown  in  figure  23. 
From  a  series  of  comparisons,  from  figures  10  and  23,  we 
found  that  initially  at  Moo=0.91  only  a  small  disturbance  is 
felt  and  expansion  waves  are  observed  at  ogive-cylinder  and 
cylinder-boattail  junctures.  Right  after  the  expansion 
waves,  the  boundary  layer  thickness  decreases  due  to  the 
increasing  speed  and  decreasing  pressure.  As  the  Mach  number 
increased  to  0.94,  shocks  are  seen  to  occur  on  the  middle  of 


(a)  Pressure  contours 
p:  0.46->0.84,  Ap=0.02 
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(b)  Mach  contours 
M:  0.85->1.0,  AM=0.01;  M:  1.0->1.25,  AM=0 .  05 


Figure  23.  Contour  plots, 


(c)  Pressure  contours 
p:  0.46->-0.86,  Ap=0.02 
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(d)  Mach  contours 
M:  0.88-1.0,  AM=0.01;  M:  1.0-1.2,  AM=0.05 


Figure  23.  (continued) 


0 


2  4 

(e)  Pressure  contours 
p:  0.46->0.86,  Ap=0.02 
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(f)  Mach  contours 
M:  0.88->1.0,  AM=0.01;  M:  1.0*1.26,  AM=0 .  02 


Figure  23.  (continued) 
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(g)  Pressure  contours 

p:  0.52-^0.92,  Ap=0.02 


(h)  Mach  contours 
M:  1.0^1.32,  AM=0.02 


Figure  23.  (continued) 


(i)  Pressure  contours 
p:  0.5-0.9,  Ap=0.02 


(J)  Mach  contours 
M:  .1.08-1.5,  AM=0.02 


Figure  23.  (continued) 
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the  cylinder  and  boattail  parts.  Right  after  the  shock 
waves,  the  boundary  layer  thickness  increases,  which  is 
clearly  shown  in  the  Mach  contour  plot,  due  to  the 
increasing  back  pressure  and  decreasing  velocity.  With  a 
further  increase  in  Mach  number,  the  regions  of  supersonic 
flow  have  grown  and  the  shocks  have  moved  downstream.  At 
Moo=0.98  speed,  the  computed  results  show  that  the  shock  in 
the  middle  of  the  cylinder  part  has  disappeared,  and  the 
only  shock  occurring  is  at  the  end  of  the  boattail  portion. 
This  is  also  evident  via  the  shear  stress  t  ►  plot  in  figure 
20.   A  further  increase  in  the  freestream  condition,  to 
M«.=1-10'  shows  the  beginning  of  a  bow  shock  at  the  nose  and 
expansion  fans  are  evident  at  the  surface  junctures.  At 
higher  speed  M<x>=1.20,  the  contour  plots  show  the  shock  as 
well  as  expansion  waves  are  more  inclined  toward  the  flow 
direction  than  the  case  of  M  =1.10. 

From  the  validity  of  the  Baldwin  and  Lomax  turbulence 
model  reported  in  reference  2,  we  believe  that  the  above 
computed  results  for  viscous  cases  are  accurate  and  good. 


CHAPTER  VIII 
RESULTS  AND  DISCUSSION:  THREE-DIMENSIONAL  FLOWS 


To  obtain  further  understanding  on  the  application  of  the 
thin- layer  Navier-Stokes  approximation,  equation  (10),  to 
three  dimensional  flow  problems,  the  projectile  at  angles  of 
attack  was  investigated. 

The  three  dimensional  grid  in  this  study  is  obtained  by 
rotating  a  hyperbolic  planar  grid  around  the  axis  of  the 
projectile  with  a  uniform  angle  increment,  A0,  in  the 
circumferential  direction.  The  angle  <f>    is  the  angle  measured 
from  the  leeward  side.  It  is  assumed  that  the  flow  problem 
considered  is  symmetric  with  respect  to  xz-plane  and  thus 
the  grid  system  generated  is  only  on  one  half  of  the 
projectile.  Two  extra  planes  are  included  to  facilitate  the 
calculation  of  the  central  difference  approximation  of  the 
spatial  derivatives  on  the  windward  and  leeward  planes. 

Based  on  experience  in  running  the  axisymmetric  code,  the 
90x60  hyperbolic  planar  grid  from  GRIDGEN  with  domain  24  is 
used  to  provide  a  good  grid  and  ensures  obtaining  accurate 
results.   This  planar  grid  was  rotated  to  form  the  three 
dimensional  grid.  The  total  grid  points  used  becomes  90,  19 
and  60  in  the  streamwise  £ ,  circumferential  n  and  radial  g 
direction,  respectively.   The  transonic  flows  of  M  =0.91. 
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0.96,  1.20  past  the  projectile  model  at  a=10°  were  chosen 
for  numerical  simulation. 


Flow  Cases  at  10°  Angle  of  Attack 
Flow  case  at  M  =0.96 

— OCT 

For  the  flow  case  Moo=0.96,  the  computed  results  are  in 
excellent  agreement  with  experimental  data  for  the  surface 
pressure  distribution.  Figures  24  shows  the  pressure 
distribution  in  the  longitudinal  direction  at  different 
circumferential  planes.  The  pressure  on  the  ogive  portion  is 
seen  to  be  higher  than  on  the  windward  side,  thus  generating 
a  upward  force  on  the  projectile.  The  boattail  region, 
however,  shows  a  reversal  with  the  higher  pressure  now 
occurring  on  the  leeward  side,  thus  causing  a  downward 
force.  The  resulting  moment  about  center  of  gravity  causes  a 
nose-up  behavior  of  the  projectile.   This  effect  can  also  be 
seen  in  the  contour  plots  of  figures  25.  The  computed 
contour  plots  also  show  that  the  shocks  on  the  leeward  side 
are  stronger  than  those  on  the  windward  side.  In  the  Mach 
contour  plot,  it  is  seen  that  the  boundary  layer  thickness 
on  the  leeward  side  is  much  thicker  than  on  the  windward 
side,  especially  after  shock  which  is  imbedded  in  the  middle 
of  boattail  part  of  leeward  side.  The  shock  positions  on  the 
windward  side  appear  to  be  slightly  downstream  of  those  on 
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Figure  24.  Computed  surface  pressure  at  various  *-statinn<= 
of  SOCBT  projectile  at  o=10°.  stations 
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Figure   24.     (continued) 
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the  leeward  side.  The  contour  plots  at  45°-135°  plane  and 
90°-90°  plane  are  shown  in  figure  26  and  27,  respectively. 
From  a  sequential  comparison,  it  is  evident  that  the 
boundary  layer  thickness  from  leeward  to  windward  becomes 
thinner  and  the  shocks  strength  also  become  weaker. 

Figures  28  shows  the  computed  pressure  distribution  and 
experimental  data  on  the  circumferential  direction  at 
different  stations.  In  these  figures,  it  is  evident  that  the 
computed  results  slightly  deviate  from  experimental  data  on 
the  boattail  portion  of  the  projectile.  These  deviations 
become  large  near  the  leeward  side  of  the  boattail  portion 
and  from  the  computed  shear  stress  t_  -distribution,  figure 
29,  it  is  seen  that  the  viscous  interaction  with  the  flow 
field  causes  reversed  flow  and  rollup  on  the  leeward  side  in 
the  vicinity  of  the  boattail.  The  reason  or  reasons  for 
these  large  deviations  are  not  known  for  certain;  potential 
sources  include  the  neglect  of  circumferential  viscous  terms 
within  the  thin-layer  viscous  model,  uncertainties  in 
three-dimensional  turbulence  modeling,  or  the  fact  that  the 
computational  projectile  model  is  different  from  that  of  the 
experiment,  since  we  extended  the  boattail  part  to  meet  a 
horizontal  sting  in  order  to  avoid  the  base  flow  region. 
Still  more  potential  sources  include  the  incapability  of  the 
two-layer  eddy  viscosity  turbulent  model  for  separated  flow 
region  which  has  been  reported  in  reference  2,  or  the  large 
truncation  error  in  the  ri-direction  A0=11.25°  or  possibly 
error  in  the  experimental  data.  This  reversed  flow  region 
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Figure  28.  Computed  surface  pressure  in  circumferential 
direction  at  different  x  station  of  SOCBT 
projectile  at  o=10° . 
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Figure   28.     (continued) 
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Figure  29.  t^  computed  at  various  ^-stations 
of  SOCBT  projectile  at  o=10° . 
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starts  to  develop  right  after  the  cylinder-boattail  juncture 
and  spreads  to  cover  the  range  of  about  -75°<0<75°  at 
midpoint  of  the  boattail  and  then  reduces  to  about 
-65°<0<65°  at  the  end  of  the  boattail. 

Flow  Cases  at  M_=0.91  and  1.20 

For  the  flow  cases  Moe=0.91  and  1.20,  the  computed  results 
are  also  accurate  as  evidenced  from  the  agreement  between 
the  computed  and  measured  surface  pressure,  as  shown  in 
figure  30.  The  distribution  of  calculated  shear  stresses, 
figure  31,  show  that  there  are  two  very  small  regions  of 
reversed  flow  in  ^-direction,  one  at  a  short  distance 
downstream  of  the  ogive-cylinder  juncture  and  the  other  at 
the  cylinder-boattail  juncture  for  1^=0.91;  however,  there 
is  none  for  the  case  of  Mao=1.20.  The  results  also  show  that 
reversed  flow  in  the  circumferential  direction  exists  on  the 
boattail  section  for  both  M  =0.91  and  1.20  and  that  the 

00 

reversed  flow  region  decreases  with  increasing  freestream 
Mach  number.  The  following  tables  list  the  relative 
importance  of  the  forces  and  moments  due  to  pressure  and 
shear  stress: 
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Figure  30.  Computed  surface  pressure  at  various  *-station< 
of  SOCBT  projectile  at  a=10° . 
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Figure    30.     (continued) 
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Figure   30.     (continued) 
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Figure   30.     (continued) 
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Figure  31.  t^  computed  at  various  ^-stations 
of  SOCBT  projectile  at  a=10° . 
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Table  3.  Forces  and  moments  due  to  pressure 
and  shear  stress. 


M 

00 

FpA 

FfA 

FpN 

FfN    <' 

<103) 

0.91 

13.21 

18.11 

126.4 

3.771 

0.96 

31.45 

20.12 

136.1 

4.142 

1.20 

89.37 

32.08 

341.8 

6.788 

N 

oo 

dp 

Df 

D 

LP 

Lf 

L    (*103) 

0.91 

34.96 

18.49 

53.45 

122.2 

0.568 

122.8 

0.96 

54.61 

20.53 

75.14 

128.6 

0.586 

129.2 

1.20 

147.4 

32.77 

180.2 

321.1 

1.115 

322.2 

H 

oo 

"p 

Mf           . 

M    (xlO3) 

0.91 

288.6 

0.747 

289.3 

0.96 

348.0 

0.780 

348.8 

1.20 

436.3 

0.733 

437.0 

The  computed  aerodynamic  forces  acting  on  the  6-caliber 
SOCBT  projectile  model  at  a=10°  show  that  the  contribution 
of  the  viscous  force  is  very  important  to  the  axial  force 
component  F^»  and  the  drag  force  component  D^ ,  since  the 
viscous  force  is  acting  along  the  surface;  however,  it  has 
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negligible  effect  on  the  normal  force  component  FfN,  the 
lift  force  component  Lf  and  the  moment  Mf  about  the  center 
of  the  projectile.  For  the  axial  force  component,  the  ratio 
of  the  pressure  part  to  the  viscous  part  F  A/FfA  is, 
respectively,  0.73,  1.53  and  2.79  for  Mw=0.91,  0.96  and 
1.20.   Several  relevant  figures  corresponding  to  Mach  number 
0.91  and  1.20  are  presented  in  figures  32  and  33  for 
reference. 


Flow  Case  at  45°  Angle  of  Attack 

For  the  application  of  the  three  dimensional  thin-layer 
Navier-Stokes  code  to  a  more  complex  flow  problem,  we  have 
considered  the  transonic  flow  of  M  =0.96  past  the  SOCBT 
projectile  model  at  a  very  high  angle  of  attack,  a=45° .  The 
same  three  dimensional  axisymmetric  grid  used  for  the  cases 
of  a=10   is  employed  and  a  converged  steady  state  solution 
has  been  obtained  with  EI=3eE=48At.   Longitudinal 
distributions  of  surface  pressure  computed  at  9 
circumferential  stations  are  presented  in  figure  34.   It 
shows  that,  in  contrast  to  the  case  of  a=10° ,  there  is  no 
shock  wave  on  the  cylinder  portion  of  the  projectile.  This 
is  also  evident  on  contour  plots  in  figure  35.  From  Mach 
contour  plot,  we  observed  that  the  local  Mach  numbers  in  the 
flow  field  around  the  body  surface  are  less  than  1,  which 
indicates  that  there  are  no  shocks  and  expansion  waves 
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Figure    32.    Computed    surface   pressure    in   circumferential 
direction   at   different   x    station   of    SOCBT 
projectile    at   a=10° . 
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Figure   32.     (continued) 


188 


-    113    - 


M  -   1.2 

ALPHA   =    10 


ii 


-  ?  -  • 


"**.—.» 


«?sstr -o ff 


.4- 


-.6 


X 


— o~. 

X 


X 


45 


JU-- 


D 


o 

D 

-x' 


30 


145 


o 

..0 

X 


186 


X 
0.89 
2.79 
3.13 
3.56 
4.22 


—  * 
...  0 

—  x 

—  □ 

—  O 


M  =    1.2 

ALPHA   =    10 


«5  30  145 

Figure   32.    (continued) 


188 


-    114   - 


M  =    1.2 

ALPHA   =    10 


Figure   32.     (continued) 
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Figure  34.  Computed  surface  pressure  at  various  ^-stations 
of  SOCBT  projectile  at  a=45° . 
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occurred  near  the  body  surface.   The  computed  shear  stress 
distribution,  x^,  as  shown  in  figure  36,  indicates  that 
there  is  no  reversed  flow  in  the  streamwise  direction; 
moreover,  it  confirms  that  there  is  no  shock  wave  on  the 
cylinder  portion.  However,  the  computed  circumferential 
shear  stress,  t,  #  distribution,  as  shown  in  figure  37, 
reveals  that  reversed  flow  exists  in  the  circumferential 
direction  over  the  entire  6-caliber  projectile  length.  As 
expected,  as  the  angle  of  attack  is  increased,  the  leeward 
side  reversed  flow  structure  becomes  stronger  and  extends 
progressively  farther  forward  along  the  projectile.  The 
appearance  of  a  reversed  flow  region  is  further  indicated  in 
figure  38,  which  shows  the  location  of  zero  shear  stress 
which  is  determined  from  interpolation  of  the  computed 
circumferential  shear  stress,  t _  ,  on  the  projectile 
surface.  This  plot  identifies  the  size  of  the  reversed  flow 
region.  The  reversed  flow  field  is  also  evidenced  in  the 
computed  results  of  the  velocity  vectors  and  streamlines. 
The  projections  of  the  velocity  field  in  a  cross-sectional 
plane  perpendicular  to  the  projectile  axis  and  its 
corresponding  pseudo  streamlines,  located  at  X/D=5.035, 
right  after  the  juncture  of  cylinder-boattail,  are  shown- in 
figure  39.   The  accuracy  of  the  results  obtained  for  this 
flow  problem  cannot  be  assessed  quantitatively,  since  there 
is  no  measurement  available  for  the  projectile  at  a=45° . 
Also,  it  has  been  reported  that  finer  grid  resolution  in  the 
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Figure  37.  x^  computed  at  various  ^-stations 
of  SOCBT  projectile  at  a=45°. 
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Figure  38.  Location  of  zero  circumferential  shear  stress 
for  SOCBT  at  a=45° . 
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Figure  39.  A  sample  velocity  field  and  trajectories  in  a 
cross  plane  at  a  short  distance  downstream  of 
the  cylinder-boattail  juncture. 
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circumferential  direction,  A$=2.50,  is  needed  for  accurate 
prediction  of  the  cross  flow  field  [5].   However,  based  on 
the  published  results  examining  the  validity  of  the 
Baldwin-Lomax  turbulence  model,  we  believe  that  the  results 
computed  and  presented  are  at  least  qualitatively  correct. 


Comparison  of  Flow  Cases  at  Different  Angles  of  Attack 

For  the  case  of  1^=0.96,  the  computed  forces  and  moment 
are  summarized  in  the  following  table  for  angles  of  attack 
of  0°,  10°  and  45°: 


Table  4.  Forces  and  moment  acting  on  the  projectile 
at  different  angles  of  attack. 


a 

°P 

Df 

D 

L 

M  (*103) 

0° 

40.28 

20.53 

60.81 

0 

0 

10° 

54.61 

20.53 

75.14 

129.2 

348.8 

45° 

1564. 

17.65 

1582. 

1608. 

761.1 

We  observed  that  there  is  a  sudden  change  of  drag  force 
between  angles  of  attack  of  10°  and  45°,  which  mainly  comes 
from  the  pressure  drag,  this  is  because  the  flow  field  on 
the  leeward  side  of  the  projectile  is  entirely  immersed  on 
the  reversed  flow  region  at  45°  angle  of  attack.  In  the 
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reversed  flow  region,  the  pressure  is  smaller  in  comparison 
with  the  pressure  on  the  windward  side.  Consequently,  the 
drag  force  of  projectile  at  45°  is  dominated  by  the  pressure 
difference  between  the  windward  and  leeward  sides.  This  also 
causes  the  rapid  change  in  the  lift  force. 


CHAPTER  IX 
CONCLUDING  REMARKS 


Numerical  simulation  of  transonic  flow  over  an  SOCBT 
projectile  configuration  for  both  the  axisymmetric  flow 
cases  and  at  angles  of  attack  have  been  obtained  using  the 
thin-layer  Navier-Stokes  codes.  Extensive,  detailed 
comparisons  of  the  computed  results  to  experimental  data  at 
Moo=0.96  have  been  made  to  evaluate  the  applicability  and 
accuracy  of  the  numerical  technique.  The  comparisons  between 
the  computed  and  experimental  surface  pressure  distributions 
show  excellent  agreement.  At  angles  of  attack,  however,  the 
solution  could  be  degraded  on  the  leeward  side  of  the 
boattail  part  due  to  the  inability  to  model  accurately  the 
developing  leeward  side  vortex  flow.  The  precise  cause  of 
the  inaccuracy  has  not  yet  been  determined;  however,  the 
turbulence  model  and  inclusion  of  the  circumferential 
viscous  terms  have  been  identified  as  areas  for  further 
investigation. 

The  excellent  agreement  between  computational  and 
experimental  results  for  the  surface  pressure  distribution 
indicated  an  accurate  computation  of  the  shear  stress 
distribution.  The  computed  results  show  that  the  viscous 
drag  is  as  important  as  the  pressure  drag  in  the  lower 
transonic  flow  regime;  and  that  it  is  less  important  when 
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Mach  number  is  increased.  At  the  same  Mach  number,  drag  due 
to  the  pressure  is  important  for  increasing  angles  of 
attack,  since  the  inclusion  of  wave  drag  and  the  pressure 
intensities  are  smaller  on  the  leeward  side  when  compared 
with  the  windward  side.   However,  the  drag  due  to  viscous 
effect  becomes  less  important.   These  results,  as  expected, 
indicate  that  the  thin-layer  Navier-Stokes  codes  employed  in 
this  study  represent  a  viable  computational  tool  for 
predicting  flow  fields  around  projectiles  at  high  angles  of 
attack. 


APPENDIX 


The  Jacobian  matrices,  A,  B  and  C,  which  appeared  in 
equation  (22)  are 


A,  B  or  C  = 


kx0   -  u9 
Kytf2  -  v9 

kz<i>2  -   we 


Kt  +  8  -  icx(y-2)u 
KXV  "  ICy(y-1)U 

Kxw  '  Kz(y-2)w 


-9(yE/P-20z)  Kx(yE/P-*^)-(y-i)u8 


Kyu  "  Kx(y"1)v 
Kt  +  6  -  k  (2T-2)v 

kvw  -  K„(y-i)v 


k2u  -  Kx(y-i)w 

KZV  -  Ky(y-l)W 

Kt  +  8  -  Kz(2f-2)w 


Ky(2TE/p-0z)-(y-l)v8   Kz(rE/P-02)-(y-l)w9 


*z<*-d 


Kt  +  re 


where 


8  =  kxu  +  k  v  +  kzw 


2     2     2 

uz  +  vz  +  wz 


^  =  (y-i)(- 


with  k   =   %,    or  ti  or  z,    for  A,  B  or  C  respectively. 
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The  viscous  flux  Jacobian  is 


M  = 


0       0  0  0       0 

m21   al6^(1/p)   a26C(1/p)   a36c(1/p)    0 


31   «26C(1/p)   a46C(Vp)   a56c(l/p)    0 
41   «36c(Vp)   *5*c(l/p)   o66c(l/p)    0 

m52         m53         m54   °06C(1/p), 


where 


21  =  ai6c("u/p)  +  a26?(-v/P)  +  o36^(-w/P) 

31  =  a26c("U/p)  +  a46c("v/p)  +  a56c(_w/p) 

41  =  a36^(_u/p)  +  a55^("V/p)  +.  a66c(-w/p) 

2 

51  =  ai6^(_u  /p)  +  a26£(~2uv/p)  +  a35^(-2uw/p) 

«45c(-v2p)  +  a66c(-w2/p)  +  a58^(-2vw/p) 

a06c(-E/P2)  +  aQ6c[(u2  +  v2  +  W2)/P] 


m52  _  "  m21  "  anMu/p) 


m54  ~  "  m41  "  an6^(w/p) 


v0uc 
'06C 


m 


53 


m 


31  -a06C<v/p) 


aQ  =  ZV    Pr_1(Cx2  +  £y2  +  Cz2) 
0l  =  y[(4/3Kx2  +  cy2  +  Cz2] 
a4  =  y[Cx2  +  (4/3  Ky2  +  Cz2] 
a6  =  VUX2  +  C  2  +  (4/3  K  2] 


«2  =  (y/3)CxCy 
«3  =  (U/3)?XCZ 
a5  =  (y/3)CyCz 
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